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PROJECT  SUMMARY 


Over  the  last  several  years  a  realrneed  has  emerged  for  numerical  schemes  to  predict 
hypersonic  flows  around  three-dimensional  configurations.  ^ftecentlyVjhe  Army  has 
shown  considerable  interest  in  applying  3-D  PNS  schemes  to  study 

supersonic/hypersonic  flowfields  around  projectiles  and  missiles.  Although  the  existing 
noniterative  PNS  schemes  offer  an  attractive  combination  of  accuracy  and  affordability, 
such  schemes  suffer  from  several  numerical  instability  and  inaccuracy  problems,  'Trytfiis  eJ>-ia 
SB1R  Phase  I  effort  we  have  developed Tt  new  3-D  PNS  scheme  with  significantly  en¬ 
hanced  stability  and  accuracy  characteristics.  This  work  was  done  by  VRA,  Inc.,  for 
US  Army  under  contract  number  DAAL03-87-C-0015.  Dr.  Clark  H.  Lewis  was  the 
principal  investigator  for  this  work,  and  Dr.  Thomas  L.  Doligalski  of  the  Army  Research 
Office,  Durham,  NC,  served  as  the  contract  monitor. 

In  this  SB1R  Phase  I  study  a  new  three-dimensional  fully-iterative  PNS  scheme  has 
been  developed  to  study  viscous  hypersonic  flows  at  large  angles  of  attack.  The  three- 
7>  fh  ,  dimensional  PNS  scheme  developed -under  this  SBIR  Phase  I  effort>is  inherently  stable 
in  the  subsonic  as  well  as  the  supersonic  flow  regions  and,  thus,  does  not  require  any 
sublayer  approximation.  'Furthermore?  it  uses  a  generalized  PNS  formulation  to  treat 
perfect-gas  and  equilibrium-air  gas  models  in  a  unified  manner.  A  second-order 
smoothing  approach  is  used  to  damp  the  solution  oscillations,  while  a  pseudo-unsteady 
approach  is  used  to  dramatically  improve  the  solution  efficiency  without  compromising 
the  solution  accuracy.  A  new  fully-implicit  and  crossflow-coupled  shock-fitting  ap¬ 
proach  has  been  developed,  along  with  a  new  predictor-corrector  solution  scheme  to 
treat  large  crossflow  separated  regions.  •'  ' '  V/  '  tf*1 ./«  •  ■'■•••*  - 

/  >  V  \ 

In  order  to  test  this  new  3D-PNS  scheme  we  studied  the  flow  around  a  15  deg 
sphere-cone  configuration  at  a  Mach  number  of  10.6  and  an  angle  of  attack  of  20  deg. 

This  test  case  corresponds  to  one  of  the  wind-tunnel  tests  conducted  by  Cleary  (1969). 

Four  different  calculations  were  done  to  evaluate  the  effects  of  grid  refinement  and  gas 
model,  and  the  numerical  predictions  were  compared  with  the  available  experimental 
data  on  wall  heat-transfer  rates.  The  results  are  very  encouraging  and  fully  substantiate 
the  accuracy,  efficiency  and  stability  claims  of  the  new  fully-iterative  three-dimensional 
PNS  scheme. 
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NOMENCLATURE 


a  =  local  speed  of  sound 

CFS  =  streamwise  skin-friction  coefficient 

CFW  =  crossflow  skin- friction  coefficient 

Cp  =  specific  heat  at  constant  pressure 

h  =  static  enthalpy  of  the  mixture 

J  =  determinant  of  the  transformation  Jacobian 

k  =  mixture  thermal  conductivity 

M  =  Mach  number 

n  =  iteration  number 

P,p  =  static  pressure 

PHI  =  circumferential  angle  measured  from  the  windward  side,  0 
PINF  =freestream  static  pressure,  p  ^ 

Pr  =  Prandtl  number 

PW  =  wall  pressure 

QW  =  total  wall  heat-transfer  rate  in  Btu/ft  2/sec 
RB  =  local  body  radius 

Re  =  Reynolds  number,  (p'V'Rn*)/^’ 

RN,Rn  =  nose  radius 

RSHK  =  radial  location  of  the  bow  shock 


T  =  static  temperature 

t  =  time 

u  =  x-component  of  mass-averaged  velocity 

Uj  =  u,  v  and  w  for  j  =  1 ,2  and  3 

Uj  =  contravariant  velocities,  uk£JiXk 

UINF  =  freestream  velocity,  V  ^ 

V  =  total  mass-averaged  velocity 

v  =  y-component  of  mass-averaged  velocity 

w  =  z-component  of  mass-averaged  velocity 

X,x  =  coordinate  along  body  axis 

Xj  =  x,  y  and  z  for  j=  1,  2  and  3 

a  =  angle  of  attack 

y  =  specific-heat  ratio 

£  *MJReM 

H  -  mixture  viscosity 

<*; 1  =  marching  or  streamwise  coordinate 

£2  =  coordinate  measured  from  the  body  to  the  outer  bow  shock 

£  j  =  coordinate  measured  from  the  windward  to  the  leeward  direction 

p  -  mixture  density 

4>  =  circumferential  angle  measured  from  the  windward  side 

Superscript 

j  =  index  in  <*;  direction 

n  =  index  for  iteration 

T  =  vector  or  matrix  transpose 

*  =  dimensional  quantity 


Subscript 

v 


L  - 


oo 


=  represents  partial  derivative 
=  freestream  quantity 
i  =  represents  i-th  chemical  species 

j,k,l  =indicial  notation  representing  1,  2  and  3 

s  =  shock  quantity 

w  =  wall  quantity 

Vector  and  Matrix  Notation 
Vector  =  bold  lower-case  character 

Matrix  =  bold  upper-case  character 

•  =  vector  dot  product 


I.  INTRODUCTION 


Over  the  past  several  years  significant  advances  have  been  made  in  the  field  of 
Computational  Fluid  Mechanics  (CFD).  Several  recent  advances  in  computer  facilities, 
their  speed  and  storage  capabilities  have  made  it  possible  to  develop  accurate  and  so¬ 
phisticated  numerical  prediction  schemes  for  various  high-speed,  fluid-flow  problems. 
On  the  other  hand  the  development  of  high-speed  experimental  facilities  and  their  oper¬ 
ation  has  remained  as  an  extremely  expensive  option. 

The  existing  methodology  for  external-flow  prediction  over  ballistic  geometries  con¬ 
sists  of  Navier-Stokes  (NS),  Parabolized  Navier-Stokes  (PNS)  and  Viscous  Shock-Layer 
(VSL)  schemes.  The  NS  schemes  (Sahu,  1986;  Sahu,  1987,  Sahu  and  Nietubicz,  1987; 
Richardson,  1986;  Kumar,  1986;  and  Rizk  et  al.,  1981)  are  typically  very  time  consuming 
and  not  well  suited  for  various  parametric  studies  required  for  design  and  analysis  pur¬ 
poses.  On  the  other  hand,  the  existing  noniterative  PNS  schemes  (Weinacht  et  al.,  1986; 
Guidos  and  Sturek,  1987;  Kaul  and  Chaussee,  1983;  and  Stalnaker,  1985)  suffer  from 
instabilities  and  inaccuracies.  Apart  from  being  noniterative,  such  PNS  schemes  require 
a  substantial  approximation  in  the  way  the  subsonic  sublayer  region  is  treated.  There 
are  generally  large  global  conservation-of-mass  (as  well  as  momentum  and  energy)  errors 
associated  with  these  PNS  methods,  which  originate  from  the  basic  solution  scheme. 
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The  existing  VSL  schemes  (Thareja  et  al.,  1983a  and  1983b;  Thompson  et  al.,  1983;  and 
Bhutta  et  al.,  1985b)  have  a  basic  limitation  of  being  parabolic  in  the  crossflow  direction 
and,  consequently,  can  not  march  through  crossflow  separated  regions.  This  prevents 
the  VSL  scheme  from  accurately  predicting  the  aerodynamic  response  of  complex 
projectile/missile  configurations,  which  may  experience  strong  crossflow  separation  ei¬ 
ther  due  a  large-angle-of-attack  condition  or  simply  due  the  3-D  nature  of  the  geometry 
being  considered  (such  as  finned  configurations,  etc.).  However,  even  under  these  con¬ 
ditions  the  flowfield  in  the  nose  region  is  attached,  and  the  VSL  schemes  represent  an 
accurate  and  efficient  way  of  generating  the  nose  solutions  for  starting  other  (more  ac¬ 
curate)  afterbody  methods  which  can  treat  crossflow  separation  (such  as  the  PNS 
schemes). 

In  this  Phase  I  effort  we  have  extended  our  basic  PNS  scheme  to  study  three- 
dimensional  supersonic/hypersonic  flowfields  using  a  general  curvilinear  coordinate  sys¬ 
tem.  However,  only  axisymmetric  geometries  have  been  considered  in  this  Phase  I  study. 
Treatment  of  complex  three-dimensional  geometries  requires  additional  considerations 
related  to  grid  generation  and  geometry  modeling,  and  will  be  addressed  in  the  follow-on 
Phase  II  effort.  In  this  Phase  I  study  substantial  effort  was  also  devoted  to  the  devel¬ 
opment  of  a  new  fully-implicit  and  crossflow-coupled  shock-fitting  scheme  in  which  the 
bow  shock  shape  is  predicted  as  the  solution  marches  down  the  body.  Furthermore,  a 
new  predictor-corrector  solution  scheme  was  developed  and  used  to  treat  the  strong 
crossflow  coupling  effects  in  and  around  the  crossflow  separated  region.  This  3-D 
perfect-gas  PNS  scheme  was  further  extended  to  include  an  equilibrium-air  gas  model. 

The  3-D  PNS  scheme  developed  under  this  Phase  I  effort  can  treat  flows  over  a  Mach 
number  range  of  3  to  30  and  for  freestream  Reynolds  numbers  as  low  as  103  (based  on 
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nose  radius).  The  present  scheme  has  been  tested  for  angles  of  attack  up  to  20  deg; 
however,  it  should  also  be  applicable  at  higher  angles  of  attack  as  long  as  the  flow  re¬ 
mains  axially  attached.  Although  this  Phase  I  study  was  restricted  to  laminar  flows,  the 
solution  scheme  can  be  easily  extended  to  include  turbulent  flow  modeling. 

This  final  report  discusses  in  detail  the  various  mathematical  and  numerical  develop¬ 
ments  done  under  this  Phase  I  effort.  One  large  angle  of  attack  test  case  was  considered 
to  validate  and  demonstrate  this  3-D  PNS  capability.  This  test  case  deals  with  the  flow 
over  a  sphere-cone  vehicle  at  a  Mach  number  of  10.6  and  an  angle  of  attack  of  20  deg. 
These  conditions  correspond  to  one  of  the  test  cases  studied  by  Cleary  (1969),  and  was 
chosen  because  in  this  case  most  of  the  flow  was  laminar  (except  for  some  turbulence 
effects  on  the  lee  side).  Both  perfect-gas  and  equilibrium-air  gas  models  were  used  to 
study  this  test  case.  Perfect-gas  flowfield  predictions  were  done  using  three  different 
computational  grids  (Cases  I,  2,  and  3),  and  the  results  were  compared  with  the  exper¬ 
imental  predictions.  The  equilibrium-air  flowfield  predictions  were  done  using  a  fine 
crossflow  grid  (Case  4),  which  was  similar  to  the  fine  grid  used  for  the  corresponding 
perfect-gas  calculations  (Case  1).  The  results  of  these  numerical  tests  are  very  encour¬ 
aging  and  demonstrate  quite  well  the  various  salient  features  of  the  new  3-D  PNS 
scheme. 
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II.  GOVERNING  EQUATIONS 


2.1.  Coordinate  System 


The  coordinate  system  used  for  the  present  3-D  PNS  scheme  is  a  general  curvilinear 
coordinate  system  shown  in  Fig.  1.  Also,  a  body-fixed  orthogonal  (Cartesian) 

coordinate  system  is  chosen  such  that  the  origin  of  the  Cartesian  coordinate  system  is 
at  the  tip  of  the  blunt  nose,  and  the  x-axis  is  aligned  with  the  axis  of  the  body.  The  z- 
axis  is  chosen  as  pointing  downward  such  that  the  windward  surface  of  the  vehicle  is  on 
the  positive  z-axis  side  (see  Fig.  1).  For  convenience  of  notation  we  will  also  refer  to  the 
x  coordinate  as  x,,  the  y  coordinate  as  x2  and  the  z  coordinate  as  x3. 

The  <*;,  coordinate  is  along  the  body  and  is  also  the  marching  direction.  The  co¬ 
ordinate  stretches  from  the  body  to  the  outer  bow  shock  and  lies  in  an  axis-normal 
plane.  The  coordinate  is  measured  in  the  crossflow  direction  from  the  windward  pitch 
plane.  In  general,  it  is  assumed  that  the  (x,y,z)  space  is  uniquely  transformable  to  the  ( 
<J3)  space  through  relations  of  the  form 
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l,  =  Zi  (X,y,z) 


(2.1) 


The  uniqueness  property  is  automatically  satisfied  as  long  as  coordinate  lines  of  the 
same  family  do  not  cross,  and  it  is  important  that  the  inverse  transform  of  Eq.  (2.1)  is 
definable.  Suppose  we  denote  the  marching  step  at  which  we  seek  the  solution  as  'j  +  1' 
and  the  previous  step  as  'j',  then  the  transformation  of  Eq.  (2.1)  is  chosen  such  that  the 
physical  (x,y,z)  grid  between  'j'  and  'j+  1'  steps  transforms  to  a  rectangle  in  the  compu¬ 
tational  (£„  £2,  £3)  plane.  The  body  surface  corresponds  to  the  <J2  =  0  curve,  whereas  the 
outer  bow  shock  corresponds  to  <S;2=LMAX  curve  (LMAX  being  the  number  of  grid 
points  in  the  <J2  coordinate  direction).  In  the  windward  pitch  plane  £3  =  0,  and  in  the 
leeward  pitch  plane  £3=KMAX,  where  KMAX  is  the  number  of  circumferential  grid 
points.  Also,  <J[  =  0  at  the  j-th  step  and  =  1  at  the  j  +  1  step.  Thus,  at  each  marching 
step,  every  grid  cell  in  the  x-y-z  space  between  the  j  and  the  j  + 1  step  is  transformed  into 
a  unit  cube  in  the  ,<J2,<*3  plane  with  =  A£2*  A£j=  1  (see  Fig.  2). 

The  transformation  given  by  Eq.  (2.1)  is  generally  difficult  to  obtain.  However,  if 
we  assume  for  the  present  that  the  (x,y,z)  grid  at  the  j  +  1  step  were  known  (subsequent 
sections  will  discuss  this  in  more  detail),  then  the  metric  derivatives  for  the  inverse 
transform 

can  be  easily  obtained  numerically.  At  each  grid  point,  this  information  about  the 
inverse-transform  metrics  is  used  to  determine  the  transformation  Jacobians  (J)  and  the 
metric  derivatives  (<?iiXj)  for  the  transformation  given  by  Eq.  (2.1). 
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2.2.  Governing  Equations 


The  full  Navier-Stokes  equations  governing  the  three-dimensional  flow  problem  can 
be  written  in  a  nondimensional  form  as  (see  Appendix  A) 


(ej  -£g,).x,  =  P  (2.3) 

We  choose  our  unknowns  to  be  the  density  ( p ),  the  density-velocity  products  (pu, 
pv  and  pw),  the  density-temperature  product  (pT)  and  the  pressure  (p).  Thus  our  vector 
of  unknowns  is 

q  =  lp,  pu,  pv,  pw,  pT,  p]T  (2.4) 

Following  the  approach  of  Peyret  and  Viviand  (1975),  it  can  be  shown  that  Eq.  (2.3) 
can  be  transformed  into  the  general  curvilinear  coordinate  system  (^),  i.e., 

(fi  "esj).{j  "  h  (2-5) 

It  is  assumed  that  in  the  freestream  the  gas  mixture  behaves  like  a  perfect  gas,  and  the 
nondimensionalization  used  in  the  above  equations  is  shown  in  Appendix  A. 


Equation  (2.5)  is  elliptic  in  £)t  Z2  and  £3  directions.  If  we  neglect  the  viscous  dissi¬ 
pation  effects  in  the  direction,  we  can  combine  Eqs.  (2.5)  and  (2.3)  in  the  following 
vectorial  equation: 


fl.tj  ~ES  l(2 


-**  3,(j  =  h 


(2.6) 
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where  the  various  components  of  this  vectorial  equations  are  defined  in  Appendix  A. 


These  five  equations  representing  the  differential  conservation  of  mass,  momentum 
and  energy  are  mathematically  closed  by  using  the  equation  of  state  for  the  particular 
gas  model  being  used.  This  equation  of  state  can  be  written  in  a  general  functional  form 
as 


HP.T.P)  =  0  (2.7) 

In  the  case  of  a  perfect-gas  model,  the  gas  is  assumed  to  be  thermally  as  well  as 
calorically  perfect.  The  equation  of  state  for  this  case  is  written  in  the  simple  algebraic 
form  as 


f(p,T,p)  =  yp  -  pT  =  0  (2.8) 

In  the  case  of  an  equilibrium-air  gas  model,  the  equation  of  state  is  written  in  the 
following  functional  form. 

Hp.T.p)  =  p  -  p(p,T)  =  0  (2.9) 

Here  the  mixture  density  is  a  function  of  pressure  and  temperature  and  is  available  in  the 
form  of  a  table. 


! 
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III.  NUMERICAL  SOLUTION  SCHEME 


This  chapter  discusses  in  detail  the  procedure  used  to  numerically  solve  the  governing 
PNS  equations. 


3.1. 


Delta-Form  of  the  Governing  Equations 


In  this  'delta  formulation'  we  solve  for  the  changes  in  the  flowfield  variables  from  one 
iteration  to  the  next.  The  approach  used  is  separately  discussed  for  perfect-gas  and 
equilibrium-air  gas  models. 


3.1.1.  Perfect-Gas  Model 

Let  us  denote  the  iteration  level  by  the  index  'n',  so  that  the  iteration  at  which  we 
seek  the  solution  is  represented  by  the  superscript  'n+  V,  and  the  previous  iteration  (the 
solution  to  which  is  known)  is  represented  by  the  superscript  'n'.  Thus,  for  the  'n+  V 
iteration  at  the  'j+  1'  marching  step  we  can  write  Eq.  (2.6)  as 
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(3.1) 


fj+l.n+l  __j+Ln+l  ,-j+I.n+l  _uj+l.n+l 
*  i.{j  _eS2.J2  ~CSJ.t3  ~h 

If  we  assume  that  the  solution  at  the  'n+  1'  level  is  close  to  the  solution  at  the  n-th  it¬ 
eration,  we  can  use  a  first-order  Taylor  series  expansion  around  the  previous  iteration 
to  write 

fj-M.n+l^fj+l.n  +  An<Aqn+l 
S  j+1.  n-H  a  s  H-l. n  +  m"  *  Aqn+1 

s^1>n+!  a4+1’n  +  M^.Aqn+1 

hi+1.n+l^hj+l.n  +  An#Aqn+l 

where 

Aqn+1  =qj+1,n+1  -qj+1>n 

and  the  matrices  A<,,  M2  and  M3  are  called  the  jacobian  matrices  (not  to  be  confused  with 
the  transformation  Jacobian  matrix).  The  forms  of  these  matrices  are  given  in  Appendix 
B.  It  should  be  noted  that  in  evaluating  the  jacobian  matrices  and  doing  the  Taylor  se¬ 
ries  expansion  around  the  n-th  iteration,  we  only  consider  the  flowfield  variables  as  the 
unknowns.  Although  the  grid  also  changes  from  one  iteration  to  the  next,  it  is  assumed 
that  these  changes  are  small  and  do  not  contribute  to  the  jacobian  matrices. 

Thus,  we  see  that  by  expanding  the  solution  around  the  n-th  iteration  and  using 
two-point  streamwise  differencing,  we  can  write  Eq.  (3.1)  as 

(Aj/Af,  -Ao)n  •  Aqn+1  -  [(A2  -eMj) .  Aqn+1]  .  -  [(A"  -eM?) .  Aqn+1], 

2  5  (3.4) 

- -  «U2  -  ~h]J+I’  "  =  g,+1’  ” 


(3.2) 


(3.3) 
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It  has  been  shown  by  Bhutta  and  Lewis  (1985a-d)  and  Bhutta  et  al.  (1985a)  that  for 
the  iterative  process  of  Eq.  (3.1),  the  simple  two-point  streamwise  differencing  is  con¬ 
servative  in  the  limit  of  convergence.  This  is  not  only  important  from  a  storage  point 
of  view,  but  it  also  gives  the  present  scheme  a  significantly  improved  capability  for 
treating  strong  compression  discontinuities. 

Equation  (3.1)  is  elliptic  in  the  £2  and  <f3  directions  so  that  for  second-order  accuracy 
we  use  central-differenced  approximations  for  all  £2  and  derivatives.  However,  the 
use  of  central-differenced  schemes  is  typically  associated  with  solution  oscillations 
(Bhutta  and  Lewis,  1985a-d;  Schiff  and  Steger,  1979;  and  Shanks  et  al.,  1979).  This 
oscillatory  behavior  becomes  more  pronounced  if  the  local  velocities  are  small,  so  that 
the  diagonal  terms  of  the  jacobian  matrices  become  relatively  small  also.  In  order  to 
damp  these  solution  oscillations,  it  is  necessary  to  add  some  additional  higher-order  dif¬ 
fusion  terms  to  Eq.  (3.1).  In  our  earlier  work  (Bhutta  and  Lewis,  1985a-d)  we  developed 
a  second-order  accurate  fully-implicit  smoothing  approach  which  is  accurate  and  simple 
to  use.  It  is  shown  in  Appendix  C  that  by  extending  the  basic  approach  of  Bhutta  nd 
Lewis  (1985a-d)  we  can  write  Eq.  (3.1)  as 

fttj,n+1  =£S&  n+1  +  «*£',• 0+1  +  hi+,,n+1  +  *.(q)(^2)2  +  nMmJ  (3-5) 
where  the  form  of  the  vectors  n  y  and  n  2  is  chosen  as 

{(Aj/A^ i  —  Aq)  »q  {  «  +  [(A2  —  eM2)  •  q  «  «  ]  «  + 

2  2  2  2  2  (3.6a) 

[(A,  -eM3).qiW2],{j)/4 

=  {(Ai/AiJ i  —  Aq)  •  q  4-  [( A2  — eM2)  •  q  {3{3],{2  + 

(3.6b) 

[(A3  -£M3).q>{3{j],{3}/4 
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It  has  been  shown  in  Appendix  C  that  (to  second-order  accuracy)  we  can  rewrite  Eq. 
(3.5)  in  terms  of  an  intermediate  solution  xi+1  as 

[fj(Xi+,)3,ij  =  £[s2(XJ+I)].{2  +  £[s3(Xi+1)]>{3  +h(xi+1)  +  0(A£2)2  +  0(AQ2  (3.7) 

The  actual  solution  that  we  seek  at  the  j  +  1  step  is  related  to  this  intermediate  solution 
by 


(x*)i+1  -  Xj+1  +  X.fjfjA^/4 

(3.8a) 

|i+1=(X*)i+,+(x‘)lf2i2A^/4 

(3.8b) 

For  this  perfect-gas  model  the  mixture  viscosity  was  obtained  using  the  Sutherland 
formula  (White,  1974),  and  the  specific-heat  ratio  was  assumed  to  be  a  constant  (1.4  for 
air).  The  mixture  Prandtl  number  was  also  assumed  fixed  (0.72  for  air),  and  the  mixture 
thermal  conductivity  was  obtained  from  the  definition  of  mixture  Prandtl  number. 

3.1.2.  Equilibrium- Air  Gas  Model 

The  numerical  scheme  used  for  equilibrium-air  flows  is  essentially  the  same  as  for  the 
case  perfect-gas  case,  except  that  the  mixture  thermodynamic  and  transport  properties 
are  provided  in  the  form  of  a  table.  The  dependant  variables  for  these  tabular  data  are 
chosen  to  be  pressure  and  temperature.  The  thermodynamic  properties  involve  the 
mixture  enthalpy,  h(p,T),  and  mixture  density,  p(p,T),  data  and  are  based  on  the  tabular 
data  of  Miner  et  al.  (1971).  The  transport  properties  involve  the  mixture  viscosity  data, 
*i(p,T),  the  mixture  thermal  conductivity  data,  k(p,T),  and  the  Prandtl  number  data, 
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Pr(p,T).  The  viscosity  and  thermal  conductivity  data  are  based  on  the  data  developed 
by  Peng  and  Pindroh  (1962).  The  Prandtl  number  data  were  developed  using  this 
viscosity  and  thermal  conductivity  data  and  the  mixture  specific-heat  data  obtained  by 
numerically  differentiating  the  aforementioned  enthalpy  data  of  Miner  et  al.  (1971).  The 
equilibrium-air  thermodynamic  and  transport  property  table  generated  for  use  in  this 
Phase  I  study  covers  the  temperature  range  of  10-15000  Kelvin  and  the  pressure  range 
of  0.0025-15.849  atmospheres.  The  pressure  range  of  the  thermodynamic  property  data 
of  Miner  et  al.  (1971)  is  from  0.000025  atmospheres  to  39810.72  atmospheres;  however, 
the  reduced  limits  of  the  current  table  are  due  to  transport  property  data  of  Peng  and 
Pindroh  (1962).  The  orginal  transport  property  data  of  Peng  and  Pindroh  (1962)  are  in 
terms  of  a  table  of  density  and  temperature.  When  these  data  were  rearranged  interms 
of  pressure  and  temperature,  it  was  observed  that  out  side  the  limits  of  0.0025-15.849 
atmospheres,  the  available  data  did  not  completely  cover  the  temperature  range  of 
10-15000  Kelvin.  Since  a  direct  extrapolation  of  these  data  in  terms  of  pressures  may 
not  be  accurate,  the  current  table  was  limited  to  this  pressure  and  temperature  range. 
It  should  be  noted  that  this  range  adequately  covers  most  of  the  flight  regime  in  which 
the  equilibrium-air  effects  may  be  important. 


3.2.  Predictor-Corrector  Solution  Scheme 


Under  large  angle  of  attack  conditions  strong  crossflow  separated  regions  may  de¬ 
velop  on  the  leeward  side.  Under  these  conditions,  solution  coupling  in  the  crossflow 
direction  is  very  important.  If  these  coupling  effects  are  not  properly  considered  during 
the  iterative  solution,  they  can  cause  severe  convergence  difficulties.  In  order  to  address 
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that  problem  of  crossflow  coupling  we  have  developed  a  new  predictor-corrector  ap¬ 
proach.  This  predictor-corrector  scheme  is  divided  into  three  different  parts;  namely,  (a) 
the  predictor  step,  (b)  the  shock  solution  and  (c)  the  corrector  step.  The  following  dis¬ 
cussion  summarizes  the  details  of  these  solution  steps. 

Using  a  two-point  streamwise  differencing  and  central-differenced  approximations  in 
the  £ 2  and  £ 3  directions,  the  final  differenced  equations  corresponding  to  Eq.  (3.4)  can 
be  written  in  the  following  block-pentadiagonal  form. 

[D.Axk_u]  +  [A.Afc/-,]  +  CB*AXk/]  +  CC.Axk/+,]  +  [E*Axk+l/] 

(3.9) 


3.2.1.  Predictor  Step 

In  the  predictor  step  we  first  neglect  the  implicit  crossflow  coupling  effects  in  favor 
of  the  body-normal  coupling  effects.  With  this  assumption  the  governing  equations  for 
the  predictor  step  become 

[A.AX*k/_,]  +  [B.Ax‘k/]  +  [C.AX*k/+1]  =  g"  (3.10) 

These  equations  are  inverted  from  the  body  to  the  shock  to  develop  a  recursive  re¬ 
lationship  between  the  solution  at  each  successive  grid  point  in  the  body-normal  direc¬ 
tion.  These  recursive  relations  have  the  form 
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AXk/  -  -Rk/*AXk/+i  +  rk/  where  k-l,..,KMAX 

/=  1,..,LMAX  —  1 


(3.11) 


3.2.2.  Crossflow-Coupled  Shock  Solution 

The  crossflow-coupled  shock-fitting  scheme  is  discussed  in  detail  in  the  following 
section.  Briefly  speaking,  using  the  recursive  relation  from  the  predictor  step  at 
<?=(LMAX-1)  location,  the  Rankine-Hugoniot  shock-crossing  equations  are  solved  to 
obtain  the  solution  at  the  shock.  This  shock  solution  is  than  used  to  solve  the  corrector 
step. 


3.2.3.  Corrector  Step 


Just  like  the  shock-point  solution,  the  solution  in  the  corrector  step  uses  the  recursive 
relations  from  the  predictor  step  to  eliminate  the  (k/-l)  contributions  in  the  difference 
molecule.  Then  using  the  fact  that  the  solution  at  the  (k/+  1)  point  is  known,  one  can 
reduce  the  Eqs.  (3.9)  to  only  a  coupled  system  in  the  crossflow  direction.  This  crossflow 
corrector  solution  can  be  written  as 


D*AXk-u  +  (B -  ARk/) •  AXl(/  +  E.Axkfl/ 


(3.12) 


=  (g n  —  A  •  r  —  C»Axic/+i 
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where  it  is  assumed  that 


AXk/-i  =  Ax  ky_1  —  -  Ry.,  •  Ax  k/  +  r 


=  —  Rk/-i  *  Ax  kj  +  r  k/-l 


(3.13) 


This  implicit  crossflow  solution  is  obtained  using  plane-of-symmetry  boundary  condi¬ 
tions  applied  in  the  windward  and  leeward  pitch  planes.  In  this  way  the  flowneld  sol¬ 
ution  is  marched  from  the  shock  to  the  body. 


3.3.  Pseudo-Unsteady  Solution  Algorithm 


It  should  be  noted  that  the  right-hand  side  of  Hq.  (3.4)  as  well  as  the  right-hand  side 
of  Eq.  (3.9)  is  the  governing  differential  equation  corresponding  to  the  fluid  mechanics 
problem  written  at  the  n-th  iteration  level  and  goes  to  zero  in  the  limit  of  conveigence. 
As  discussed  by  Bhutta  and  Lewis  (1985a-d),  under  these  conditions  the  exact  form  of 
the  left-hand  implicit  terms  is  of  no  great  consequence  except  that  it  affects  the  conver¬ 
gence  path  of  the  solution.  With  this  idea  in  mind  we  do  not  update  the  Jacobian  ma¬ 
trices  beyond  the  first  iteration;  i.e.,  we  assume  that 


Ao  *  Ag-! 

A?  “AT* 

MS  *  Mf’ 
MS  a  Mr1 


(3.14) 
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Such  a  pseudo-unsteady  approach  has  been  discussed  in  depth  by  Bhutta  and  Lewis 
(1985a-d)  and  interested  readers  are  referred  to  these  references  for  further  details.  Using 
this  pseudo  unsteady  approach  we  obtain  the  final  differenced  equations  as 

-  Ac)"*1  .  AXn+i  +  [(A2  -eMj)"*1  •  Aln+%  +  C(A3  -£M3)n=1 .  AXn+,]>{j 

(3.15) 

-  “  -£Sit2  ~eht3  -h]i+1,n  =  gi+1,n 

The  converged  limit  of  Eq.  (3.15)  is  the  same  as  the  converged  limit  of  Eq.  (3.4).  How¬ 
ever,  Eq.  (3.15)  is  considerably  more  efficient  and  faster  to  solve.  With  the  present 
pseudo-unsteady  approach  the  time  for  each  iteration  after  the  first  iteration  (n=  2,3,...) 
is  only  10-15%  of  the  time  taken  by  the  first  iteration. 


3.4.  Boundary  Conditions 


The  problem  represented  by  the  governing  PNS  equations  is  a  split-boundary-value 
problem;  i.e.,  the  equations  are  hyperbolic-parabolic  in  the  £ ,  direction  and  elliptic  in  the 
£2  and  £3  directions.  Thus,  in  order  to  solve  the  problem  completely  we  need  initial 
conditions  to  be  specified  at  the  start  of  the  marching  procedure,  boundary  conditions 
to  be  specified  at  the  wall  and  at  the  outer  bow  shock  and  boundary  conditions  to  be 
specified  in  the  windward  and  leeward  pitch  planes  (for  flows  with  a  pitch-plane  of 
symmetry). 
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3.4.1.  Initial  Conditions 


The  initial  conditions  to  start  the  perfect-gas  and  equilibrium-air  PNS  solutions  were 
obtained  from  a  VSL  blunt-body  solution  scheme  for  perfect-gas  and  equilibrium-air 
flows  (Murray  and  Lewis,  1978;  Thareja  et  al.,  1983a;  and  Thompson  et  al.,  1983)  The 
quality  of  such  VSL  solutions  has  been  discussed  in  great  detail  by  Thompson  et  al. 
(1983)  and  Bhutta  et  al.  (1985b).  The  VSL  blunt-body  solution  is  interpolated  to  obtain 
the  starting  solution  at  the  initial  data  plane  (I DP)  for  the  3-D  PNS  afterbody  solution. 
We  typically  choose  the  starting  location  to  be  approximately  2-3  nose-radii  downstream 
of  the  nose  stagnation-point  location. 

3.3.2.  Wall  Boundary  Conditions 

The  boundary  conditions  at  the  wall  consist  of  six  independent  relations  representing 
the  nature  of  the  gas  mixture  and  the  physical  conditions  at  the  wall.  These  conditions 
are: 

(1)  Equation  of  state  of  the  gas:  f(p,p,T)  =  0 

(2)  No-slip  condition  for  'u'  velocity  component:  pu  =  0 

(3)  No  slip  condition  for  V  velocity  component:  pv  =  0 

(4)  No  slip  condition  for  'w'  velocity  component:  pw  =  0 

(5)  Specified  wall  temperature:  (pT)  =  pTw. 
and 

(6)  Zero  pressure  derivative  in  the  £2  direction  (p{2  =  0) 
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The  first  five  boundary  conditions  are  easy  to  visualize  as  they  represent  the  actual 
physical  conditions  at  the  wall.  The  sixth  boundary  condition  on  the  pressure  derivative 
comes  from  a  boundary-layer-type  analysis  performed  at  the  wall.  The  above  set  of 
boundary  conditions  are  well-posed  and  form  a  linearly  independent  set. 

3.3.3.  Boundary  Conditions  at  the  Shock 

The  boundary  conditions  at  the  outer  bow  shock  are,  however,  much  more  involved. 
This  boundary  condition  involves  a  fully-implict  and  crossflow  coupled  shock-fitting 
approach,  and  the  bow  shock  is  predicted  as  the  solution  marches  down  the  body.  This 
fully-implict  bow  shock-fitting  scheme  will  be  discussed  in  detail  in  Chapter  4. 

3.3.4.  Circumferential  Boundary  Conditions 

Presently,  the  proposed  three-dimensional  nonequilibrium  PNS  schemes  can  only 
treat  flows  with  a  pitch  plane  of  symmetry;  i.e,  the  vehicle  geometry  is  symmetric  with 
respect  to  the  pitch  plane  and  there  is  no  yaw.  For  such  a  case  the  boundary  conditions 
at  the  windward  and  leeward  pitch  planes  consist  of  reflective  or  symmetric  boundary 
conditions.  The  symmetric  and  reflective  boundary  conditions  used  in  the  present  study 
are  based  on  the  second-order  crossflow  boundary  conditions  used  by  Kaul  and 
Chaussee  (1983)  and  Shanks  et  al.  (1979). 
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IV.  FULLY-IMPLICIT  SHOCK-FITTING  PROCEDURE 


4.1.  Development  of  the  Shock-Fitting  Equations 


In  this  study  we  have  developed  a  new  fully-implicit  and  crossflow  coupled  shock¬ 
fitting  scheme.  In  this  scheme  the  bow  shock  location  is  iteratively  predicted  as  the 
solution  marches  down  the  body.  This  chapter  provides  a  brief  but  complete  description 
of  this  shock-fitting  approach.  The  important  features  of  this  bow  shock-fitting  ap¬ 
proach  are: 

(a)  The  bow-shock  shape  location  is  predicted  along  with  the  flowfield  solution  and  does 
not  have  to  be  specified  a  priori. 

(b)  Unlike  earlier  noniterative  shock-propagation  approaches  (Kaul  and  Chaussee,  1983; 
Shanks  et  al.,  1979;  and  Chaussee  et  al.,  1982),  the  present  shock-fitting  scheme  is 
fully  iterative  and  treats  the  various  gas  models  (perfect-gas  or  equilibrium-air  or 
nonequilibrium-air)  accurately  and  in  a  unified  manner. 

(c)  Unlike  earlier  noniterative  shock-propagation  approaches  (Kaul  and  Chaussee,  1983; 

Shanks  et  al.,  1979;  and  Chaussee  et  al.,  1982),  the  present  approach  does  not  as- 
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sume  the  flowfield  behind  the  shock  to  be  inviscid.  This  can  be  quite  important 
when  strong  flowfield  gradients  exist  behind  the  shock,  as  may  be  the  case  in  the 
nose-dominated  region  and  in  regions  where  the  bow  shock  starts  to  interact  with 
the  imbedded  shock  waves  (or  compression  waves)  originating  from  the  body. 

(d)  Unlike  the  iterative  shock-fitting  approaches  of  Helliwell  et  al.  (1980)  and  Lubard 
and  Helliwell  (1973),  the  present  shock  fitting  approach  is  not  only  for  a  general 
curvilinear  coordinate  system  but  also  does  not  increase  the  matrix  size  of  the 
block-tridiagonal  solution  between  the  body  and  the  shock. 

(e)  Unlike  any  earlier  iterative  or  noniterative  bow  shock-fitting  scheme,  the  present 
shock-fitting  scheme  does  not  neglect  the  crossflow  coupling  effects  at  the  shock. 
This  results  in  accurate  and  smooth  shock  shapes  even  when  there  are  strong 
crossflow  variations  of  the  conditions  behind  the  shock.  This  can  be  especially  im¬ 
portant  when  dealing  with  complex  3-D  configurations  (which  is  one  of  the  main 
objectives  of  this  study)  where  the  3-D  nature  of  the  body  can  interact  with  the  bow 
shock  and  substantially  distort  it.  Similar  strong  crossflow  variations  may  also  oc¬ 
cur  on  simple  configurations  when  pitched  at  large  angles  of  attack. 

In  developing  the  present  bow  shock-fitting  scheme  we  assume  that  from  one  iter¬ 
ation  to  the  next  the  shock-points  move  along  the  <j;2  grid  line.  This  direction  corre¬ 
sponds  to  the  intersection  of  the  (inconstant  and  £2  “constant  surfaces.  This 
assumption  allows  us  to  reduce  the  size  of  unknowns  to  be  solved,  and  the  final  solution 
has  only  one  additional  unknown  at  the  shock  which  completely  defines  the  spatial 
movement  of  the  shock  point.  This  smaller  size  of  the  unknowns  is  very  important  for 
a  faster  iterative  solution  and  faster  convergence  characteristics  of  the  overall  implicit 
shock-prediction  scheme.  Furthermore,  this  simplification  only  represents  a  certain 
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constraint  on  the  direction  in  which  the  shock-point  moves  and  has  no  affect  on  the 
accuracy  of  the  shock-crossing  equations. 


Furthermore,  in  order  to  simplify  the  numerical  solution,  it  is  assumed  that  the  metric 
derivatives  y  >{2  and  z  at  the  shock  can  be  safely  assumed  as  known  from  the  previous 
iteration;  i.e.. 


(y.j2).n+,-(y.{2)" 

(z  ^+1  -  (Z,/s 


(4.1) 


This  is  also  only  a  simplification  and  does  not  affect  the  accuracy  of  the  final  converged 
solution.  For  the  case  of  nonequilibrium-air  flows  a  frozen  shock-crossing  approxi¬ 
mation  is  used. 

Using  the  first  of  the  aforementioned  assumptions  we  denote  the  amount  by  which 
the  shock  point  moves  in  the  <*2  direction  as  A In  other  words 

+  (A)sn+i  (4.2) 

Thus  the  corresponding  movement  of  the  shock-point  coordinates  from  one  iteration  to 
the  next  can  be  written  as 

+  ow"(A)"+1  <4-3) 


Furthermore,  we  define  a  set  of  six  unknowns  (a  ()  such  that 
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-**+•  _ /v  \iW-l 

«1  =  (x{,)s 

«r* = (Zf.r1 

«4n+1  =  (y,{3r 


(4.4) 


Where  it  should  be  noted  that  these  five  unknowns  are  actually  functions  of  a  single 
unknown  (A,)n+1  such  that 

*1  =  *1  +  (X,{2)s“s 

*?'  -  «2  +  (y.{2)Xn+1 

«3+1  =  +  (z.{2)X+1  (4.5) 

«r’  -  «•  +  [(y.{j)X+I],{j 

«r 1  - «“  +  [(z.  i2)X+,].{3 


At  this  point  we  define  three  orthogonal  vector  directions  n,  t  and  s  at  the  shock 
surface  (note  that  these  are  not  unit  vectors).  The  vector  n  is  normal  to  the  shock  sur¬ 
face.  The  vector  t  is  tangent  to  the  shock  surface  and  directed  in  the  <J3  direction.  The 
vector  s  is  also  tangent  to  the  shock  surface  and  orthogonal  to  vectors  n  and  t 

It  can  be  shown  that  these  vectors  can  be  written  as 

2  2 

s  =  ai(a4  4-  a5)i  +  a5(a2a5  —  a3a4)j  —  a5(a2a5  —  a3a4)k 
n  =  (a2a5  -  a3a4)i  —  a,asj  +  a,a4k  (4.6) 

t  =  a4i  +  a5k 


Or 
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s  =  [sx(a„a4,a5)]i  +  [sy(a2,  a3,  a4,  ot5)]j  +  [sz(oc2,  a3,  a4,  a5)]k 
n  =  [nx(a2,  a3,  a4,  a5)]i  +  [ny^as)]!  +  [nz(a„  ot4)]k  (4.7) 

t  =  [ty(a4)]j  +  [tz(a5)]k 

With  such  a  definition  of  the  shock-normal  and  shock-tangent  vectors  the  velocity 
components  normal  and  tangent  to  the  shock  surface  can  now  be  written  as 

V,  =  (ussx  +  vsSy  +  wss2)/|s| 

Vn  =  (usnx  +  vsny  +  wsnz)/|n|  (4.8) 

Vt  =  (u,tx  +  vsty  +  wstz)/|t| 

where  V  nt  V ,  and  V ,  are  the  velocity  components  in  the  n,  t  and  s  directions,  respec¬ 
tively.  Similarly  we  can  also  define  the  freestream  velocity  components  (V  „)  (V ,)  ^ 
and  (V,)M. 

Having  defined  the  relevant  velocity  components  for  the  purpose  of  writing  the  five 
Rankine-Hugoniot  shock-crossing  equations  (representing  the  conservation  of  mass, 
momentum  and  energy),  we  note  that  we  have  actually  seven  unknowns  at  the  shock. 
These  seven  unknowns  are  written  in  a  vectorial  form  as 


qs ,  [>.  pu,  pv,  pw,  pT,  p,  A]J  (4.9) 

Thus  we  need  two  more  equations  to  close  the  system  of  equations  at  the  shock.  One 
of  these  additional  equations  is  the  equation  of  state  of  the  gas  and  the  other  equation 
is  provided  by  applying  the  differential  continuity  of  mass  equation  behind  the  shock. 
As  we  see  no  approximation  other  than  the  assumption  of  a  Rankine-Hugoniot  shock 
has  been  made.  These  equations  are  equally  valid  whether  the  conditions  behind  the 
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shock  are  viscous  or  inviscid  dominated  or  whether  substantial  flowfield  gradients  exist 
behind  the  shock. 


The  seven  governing  equations  at  the  shock  can  now  be  written  as 


h  +  V2/2  —  (ho)^  „  o 

(4.10a) 

pVn-pjyj"  =  0 

(4.10b) 

p[V«  -  (VJJ  =  0 

(4.10c) 

p[v,-(vsu  -  o 

(4.10d) 

'oo  +P(V„)2  -  Poo(Vnt  =  0 

(4.10e) 

along  with  the  equation  of  state  written  in  the  functional  form 


fl>,T,p,C,)  =  0  (4.100 

and  the  differential  continuity  equation  written  as 

(pUj/J)  j.  =  0  (4.  lOg) 

Equations  (4.10a)  to  (4.100  are  quite  straight  forward;  however,  the  differential  conti¬ 
nuity  equation  (Eq.  (4.10g)J  needs  some  further  elaboration.  Noting  the  definition  of 
contravariant  velocity  components  (U  J,  we  can  re-write  this  equation  as 

Cai,xj(Puj)],ei  =  o  (4.i0g) 


where 
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ai»  -  (z.j2a4  -  y,(2*s) 


a!y  ”  X{2«S 
a!z  =  X.i2a4 

a2x  =  (a2a5  -  a3a4) 

a2y  =  -a,a5  (4.11) 

a2z =  ala4 

a3x  =  (y,{2«3  ~  Z.{2a2) 
a3y  ==  (z, i2a\  ~  x,{2a 3) 
a3z=  (X{2a2-y.{2«l) 

It  should  be  noted  that  the  dependance  of  all  the  quantities  appearing  in  Eqs.  (4. 10a) 
to  (4.10g)  on  the  seven  unknowns  at  the  shock  point  (q,)  is  now  completely  described. 
Thus  we  can  write  the  seven  equations  (Eqs.  (4. 10a)  thru  (4. 10g)J  at  the  shock  in  the 
vectorial  form 


-  fs(q"+1)  -  0  (4.12) 

These  equations  can  now  be  linearized  around  the  previous  iteration.  Using  central- 
differenced  approximations  for  £3  derivatives  and  backward-differenced  approximations 
for  £  1  and  <J2  derivatives,  we  can  rewrite  these  equations  in  the  form 


(AS)U  •  (Aqs)k-i  +  +  (Cs)".(Aq,)kn:i 

(4.13) 


0 

tA 

O 

_ 1 

AqJo.MAX-1 1 

1 

O 

O 

- J 

1 

0 

I _ 
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I 


Here  it  should  be  noted  that  A  „  BJ  and  C ,  are  7x7  matrices,  g’  is  a  vector  of  length  7 
and  D ,  is  a  6x6  matrix.  This  equation  represents  the  final  set  of  equations  to  be  solved 
at  the  shock  point  in  conjunction  with  fiowfield  solution  within  the  shock-layer. 


4.2.  Coupling  of  the  Shock-Point  Solution  With  the  Inner 
Fiowfield 


As  can  be  seen  with  from  Eq.  (4. 1 3)  the  solution  of  the  equations  at  the  shock  are 
coupled  to  the  inner  fiowfield  solution  through  A^max-i-  1°  order  to  finally  solve  this 
system  of  equations  we  note  that  Eq.  (3.14)  shows  that  the  inner  fiowfield  solution  (from 
/=  1,2,3..,LMAX-1)  is  decoupled  in  the  crossflow  direction.  Thus  a  forward  substitution 
approach  is  used  for  all  crossflow  planes  to  develop  the  recursive  relations 

AXk/  =  +  rk/  (4-14) 

Now  we  note  that  based  on  our  smoothing  approach  we  can  write  to  a  second-order 
accuracy 

A^Xmax-i  —  AXk,LMAX-i  (4.15a) 

AR^Xm  ax  —  Ax^Xmax  (4- •  ^a) 

The  inner  fiowfield  solution  can  now  be  related  to  the  shock-point  solution  vector 
AqJ+I  through  relations  of  the  form 

Aq^LMAX— 1  =  ~  t^k.LMAX-l  ^3  *  (Aqs)k+'  +  rk,LMAX-l 
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After  substituting  Eq.  (4.16)  into  Eq.  (4.13),  we  can  reduce  Eq.  (4.13)  to  the  form 

A;-(Aq.)”!  +  BJ.(A,S);+I  +  CjJ.fAqjKj  =  gl  (4.17) 

Equation  (4.17)  is  now  solved  using  appropriate  reflective  and  symmetric  boundary 
conditions  in  the  leeward  and  windward  pitch  planes  of  symmetry.  This  solution  gives 
simultaneously  the  Aq,n+1  vectors  at  each  shock  point  (k  =  1,2,3,...,KMAX). 

Using  the  shock  point  solution  and  the  recursive  relations  of  Eq.  (4.14),  we  can  ob¬ 
tain  the  intermediate  solution  vector  Ax  n+1  for  all  interior  points.  The  final  smoothed 
solution  at  these  interior  points  is  than  determined  from  Eqs.  (3.8a)  and  (3.8b).  The  x, 
y  and  z  coordinates  of  the  new  shock  point  locations  are  now  determined  using  Eq.  (4.3) 
and  the  grid  is  updated  for  the  next  iteration.  This  overall  iterative  process  is  repeated 
until  the  solution  converges  at  all  grid  points,  and  then  the  solution  moves  on  to  the  next 
marching  step. 
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V.  A  SIMPLIFIED  ANALYSIS  FOR  DEPARTURE 


Unlike  classical  PNS  treatments  in  the  present  approach  the  governing  PNS 
equations  are  written  in  as  a  difTerential/algebraic  system  rather  than  a  pure  differential 
system.  By  a  differential/algebraic  system  we  mean  a  system  of  equations  composed  of 
several  differential  equations  coupled  together  through  an  algebraic  relation.  In  the 
present  PNS  formulation  the  differential  equations  governing  the  conservation  of  mass, 
momentum  and  energy  are  coupled  together  through  the  algebraic  relation  representing 
the  equation  of  state.  By  studying  a  model  differential/algebraic  problem  (see  Appendix 
D)  we  have  been  able  to  show  that  unlike  classical  PNS  treatments  (Schiffand  Steger, 
1979;  Shanks  et  al.,  1979;  Chaussee  et  al.,  1981;  and  Vigneron  et  al.,  1978)  the  present 
formulation  is  unconditionally  time-like  (Bhutta  and  Lewis,  1985c  and  1985d).  In  order 
to  further  demonstrate  this  approach  we  also  studied  a  simplified  version  of  the  govern¬ 
ing  PNS  equations.  The  general  outline  and  conclusions  of  this  simplified  eigenvalue 
analysis  are  briefly  presented  in  the  following  discussion. 

In  order  to  simplify  the  required  mathematics,  let  us  restrict  ourselves  to 

(a)  two-dimensional  perfect-gas  flows,  and 

(b)  an  evenly  spaced  square  grid  such  that  £u  =  £2>z  =  1  and  £Uz  =  £2jX  =  0 
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The  simplification  of  a  perfect  gas  is  done  with  the  understanding  that  even  for  the 
actual  multi  component  reacting-gas  case  the  fluid  is  still  assumed  to  be  a  mixture  of 
perfect  gases.  Consequently  the  basic  numerical  scheme  has  to  be  at  least  stable  for  a 
single-species  perfect-gas  model.  Furthermore,  we  choose  to  approximate  the  equation 
of  state  for  a  perfect  gas  by 


yp  -  pT  +  0(p  j  ,  4-  pj,3)  =  0  (5.1) 

where  the  coefficient  '0'  is  chosen  such  that  0^0,  and  for  all  practical  purposes 

6p  ti  +  0p jr  +yp  -  pT  as  yp  -  pi 


It  should  be  noted  that  the  use  of  this  coefficient  0  is  for  the  sole  purpose  of  the  fol¬ 
lowing  stability  analysis,  and  not  for  the  actual  solution  scheme.  In  other  words,  the 
actual  solution  corresponds  to  the  use  of  0  =  0.  The  reason  for  introducing  this  coeffi¬ 
cient  '0'  is  that  it  makes  the  streamwise  jacobian  matrix  nonsingular,  so  that  one  can 
perform  certain  matrix  operations  to  simplify  the  stability  analysis.  Furthermore,  there 
are  no  mathematical  tools  available  to  directly  determine  the  character  of 
differential/algebraic  systems.  The  available  mathematical  tools  are  strictly  for  purely 
differential  systems.  Thus,  by  introducing  the  coefficient  '0'  in  Eq.  (5.1)  we  are  able  to 
transform  our  differential/algebraic  problem  into  a  purely  differential  form,  which  can 
then  be  analyzed. 


The  choice  of 'small  enough  0'  does  not  adversly  affect  the  resulting  conclusions.  In 
other  words  we  can  take  the  limit  0-»O.  This  is  indeed  a  heuristic  approach,  and  the 
justification  comes  from  the  model  differential/algebraic  problem  discussed  in  Appendix 
D.  The  use  of  '0'  in  Eq.  (5.1)  corresponds  to  the  use  of  the  V  term  in  Formulation  II 
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of  the  model  problem.  For  this  model  problem  it  has  been  shown  in  Appendix  D  that 
the  choice  of  e-+0  +  does  not  produce  any  singularity  in  the  final  solution.  This  fact  has 
also  been  numerically  demonstrated  by  the  present  authors  in  their  earlier  work  (Bhutta 
and  Lewis,  1985d). 

With  the  equation  of  state  given  by  Eq.  (5.1),  and  after  neglecting  the  viscous  terms 
containing  the  contributions  of  w  and  w  j2,  we  can  write  the  simplified  PNS  equations 
as 


(A"  •  d)>{i  +  (A"  •  d)>{2  -  £[Mn  •  d]>{2  -  AS .  d 


(5.2) 


where 


d  =  Aqn+* 


(5.3) 


and  Af,  A$  and  Mn  are  the  jacobian  matrices. 

If  we  assume  that  A?,  A§  and  Mn  do  not  change  with  £,  and  (a  frozen  coefiicient 
analysis),  we  can  write 

Aj  •  d  +  Aj  •  d  j2  —  Bf  •  d  +c  =  0  (5.4) 

Although  the  above  equation  is  a  significantly  simplified  version  of  the  original  PNS 
equations,  it  is  still  difficult  to  study  directly.  As  a  further  simplification,  we  choose  to 
separately  look  at  the  viscous  and  inviscid  limits  of  Eq.  (5.4).  We  also  note  that  as  a 
minimum  criterion  of  streamwise  stability  (which  in  this  case,  and  for  the  class  of  PNS 
schemes  based  on  the  Schiff-Steger  formulation,  implies  a  marching-like  character  of  the 
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governing  equations),  the  "simplified  PNS  equations"  being  studied  have  to  be 
streamwise  stable  (i.e.,  marching-like  in  the  streamwise  direction)  in  the  viscous  as  well 
as  the  in  viscid  limits.  In  the  following  sections  we  look  into  the  streamwise  stability 
(marching-like  character)  of  the  viscous  and  the  inviscid  limits  of  Eq.  (5.4). 

5.1.  Inviscid  Limit 

The  inviscid  limit  of  Eq.  (5.4)  can  be  written  as 

«*  +  CAta-d*  +  Af,.c-0  (5.5) 

or 

d{l  +  N,.di{j  +  c  =  0  (5.6) 

This  equation  is  now  in  a  form  which  can  be  easily  studied.  For  Eq.  (5.6)  to  be  sta¬ 
ble,  the  direction  should  be  a  valid  marching  direction.  In  other  words,  Eq.  (5.6)  has 
to  be  hyperbolic/parabolic  or  marching-like  always.  This  condition  is  satisfied  if  the 
eigenvalues  of  N ,  are  all  real.  If  for  simplicity  we  assume  that  w  <  <  u,  then  an 
eigenvalue  analysis  gives  the  eigenvalues  of  N  ,  as  (Bhutta  and  Lewis,  I985d) 

=  ( l ,  w/u,  w/u,  w/u,  w/u)  (5.7) 

Thus,  we  see  that  all  the  eigenvalues  of  N  |  are  unconditionally  real.  That  is  to  say,  the 
simplified  PNS  equations  being  studied  are  unconditionally  marching-like  in  the  inviscid 
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limit  and  represent  a  stable  marching  scheme  in  the  subsonic  as  well  as  the  supersonic 
flow  regions. 

It  is  of  interest  to  consider  the  case  where  (like  conventional  non-iterative  PNS 
schemes)  we  do  not  uncouple  the  pressure  terms.  In  such  a  case  the  form  of  A  , !  is 
similar  to  the  one  studied  by  Schiff  and  Steger  (1979).  The  only  difference  is  that  we 
have  pT  as  the  independent  variable  rather  than  'e'.  The  A  , f  matrix  for  such  a  case  is 
nonsingular  and,  thus,  can  be  inverted.  However,  the  eigenvalues  of  N  ,  for  such  a  case 
(where  we  do  not  uncouple  the  pressure)  turn  out  to  be  complex  in  the  subsonic  flow 
regions,  and  the  marching  scheme  becomes  unstable  (not  marching-like;  i.e.,  elliptic) 
unless  methods  such  as  the  sublayer  approximations  of  Schiff  and  Steger  (1979)  or 
Vigneron  et  al.  (1978)  are  used. 

5.2.  Viscous  Limit 

In  the  viscous  limit,  Eq.  (5.4)  simplifies  to 

dj^tAfS^.d^  +  Af'.c  (5.8) 

or 

d.J,=N  2-d,{2?2  +  c  (5.9) 

In  this  form,  the  stability  analysis  becomes  much  simpler.  The  criterion  of  a  stable 
marching  scheme  requires  that  Eq.  (5.9)  should  be  parabolic.  The  parabolic  character 
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depends  upon  the  eigenvalues  of  N  2,  which  should  be  real.  Furthermore,  in  order  to 
have  positive  diffusion  effects  in  the  direction,  these  eigenvalues  should  also  be  posi¬ 
tive.  Thus,  for  the  viscous  limiting  case  to  be  stable,  the  eigenvalues  of  N  2  should  be 
real  and  positive.  An  eigenvalue  analysis  of  N  2  shows  that  the  eigenvalues  are  (Bhutta 
and  Lewis,  1985d) 


<Tj  =  (0, 0,  tjti/Prpu,  4ep/3pu,  ep/pu)  (5. 10) 

Thus,  the  eigenvalues  a  t  are  always  real;  however,  they  are  positive  only  if  u  >  0.  That 
is  to  say,  as  long  as  no  flow  reversal  occurs  in  the  streamwise  direction,  the  viscous  limit 
of  the  simplified  PNS  equations  is  also  unconditionally  marching-like.  Since  flow  re¬ 
versal  means  axial  separation,  this  streamwise  stability  requirement  actually  tells  us  that 
a  "single-sweep"  solution  of  these  PNS  equations  can  not  be  marched  through  regions 
of  axial  flow  separation.  Of  course,  this  conclusion  comes  as  no  surprise  and  has  been 
a  well  accepted  fact  in  fluid  mechanics  for  a  long  time.  It  may  be  of  value  to  note  that 
the  viscous  terms  do  not  include  any  pressure  terms  and,  thus,  for  the  present  scheme 
as  well  as  the  previous  PNS  schemes  the  viscous  terms  do  not  provide  any  speed-of- 
sound  contribution  to  the  eigenvalues  of  the  viscous  limit.  The  speed-of-sound  con¬ 
tributions  to  these  eigenvalues  for  the  viscous  limit  can  come  only  from  the  jacobian 
matrix  corresponding  to  the  streamwise  convective  terms.  For  the  classical  PNS 
schemes  (Schiff  and  Steger,  1979;  Shanks  et  al.,  1979;  Chaussee  et  al.,  1981;  and 
Vigneron  et  al.,  1978),  this  speed-of-sound  contribution  is  the  one  which  causes  the 
problem  of  negative  eigenvalues  in  the  subsonic  sublayer  region.  For  the  present  PNS 
scheme,  although  the  speed  of  sound  does  appear  in  the  streamwise  jacobian  matrix  ( 
Af),  it  does  not  contribute  to  the  eigenvalues  in  the  viscous  limit. 
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VI.  RESULTS  AND  DISCUSSION 


In  order  to  test  the  accuracy  and  efficiency  of  the  proposed  three-dimensional  PNS 
scheme,  we  studied  the  flow  around  a  sphere-cone  configuration  at  an  angle  of  attack 
of  20  deg  and  at  a  Mach  number  of  10.6.  This  test  case  has  been  studied  by  Cleary 
(1969)  in  the  NASA  Ames  hypersonic  wind  tunnel.  The  freestream  conditions  for  this 
case  are  given  in  Table  1  and  the  vehicle  geometry  is  shown  in  Fig.  3.  This  sphere-cone 
geometry  consists  of  a  0.028  m  nose  radius  with  a  15  deg  aftcone.  The  vehicle  geometry 
is  17.5  nose  radii  long.  In  this  study  both  perfect-gas  and  equilibrium-air  gas  models 
were  used  to  model  the  gas  chemistry. 

Cases  1,  2  and  3  represent  the  three  different  computational  grids  used  to  predict  the 
flowfield  using  a  perfect-gas  model.  Case  1  uses  31  circumferential  planes  and  50  points 
between  the  body  and  the  shock.  Case  2  uses  25  circumferential  planes  and  Case  3  uses 
17  circumferential  planes.  For  all  of  these  three  perfect-gas  calculations  (Case  1,  2  and 
3)  the  grid  generated  was  a  simple  noncylindrical  grid  based  on  equally  spaced  grid 
at  the  wall  and  at  the  shock.  These  equally  spaced  body  and  shock  points  were  joined 
using  straight  lines  in  the  £2  direction.  For  all  of  these  grids  the  £2  grid  distribution  was 
the  same  and  the  grid  spacing  near  the  wall  was  0.001%  of  the  local  shock  standoff  dis- 
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tance.  The  purpose  of  using  a  noncylindrical  grid  was  to  demonstrate  the  general- 
curvilinear-coordinate  formulation  used  in  our  3-D  PNS  scheme. 

Case  4  represents  the  corresponding  solution  using  an  equilibrium-air  gas  model. 
Like  the  Case  1  grid,  Case  4  uses  31  circumferential  planes  and  50  points  between  the 
body  and  the  shock.  However,  in  this  case  we  used  a  cylindrical  coordinate  system  in 
which  the  £3  (crossflow)  grid  points  were  located  using  equal  angular  spacing.  The  grid 
distribution  in  the  <J2  direction  was  the  same  as  the  one  used  for  the  perfect-gas  calcu¬ 
lations  (Cases  1,  2  and  3),  and  the  <*2  grid  spacing  at  the  wall  was  0.001%  of  the  local 
shock  standoff  distance. 

The  following  sections  summarize  some  of  the  important  results  of  these  studies  and 
discuss  the  solution  accuracy,  efficiency  and  grid-refinement  characteristics. 


6.1.  Perfect-Gas  Calculations 


The  following  sections  discuss  the  results  of  the  Case  1,  Case  2  and  Case  3  calcu¬ 
lations  using  a  perfect-gas  model. 

I 

6.1.1.  Predictions  of  Surface-Measurable  Quantities 

| 

Figures  4  thru  9  show  the  axial  and  crossflow  distributions  of  the  surface-measurable 
quantities  (such  as  wall  pressures,  skin-friction  coefficient  and  wall  heat-transfer  rates) 
■  using  a  perfect-gas  model  and  the  fine  crossflow  grid  (Case  1).  The  axial  distribution 
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of  wall  pressures  is  shown  in  Fig.  4,  while  the  corresponding  crossflow  distribution  is 
shown  in  Fig.  5.  These  results  show  that  the  streamwise  and  crossflow  distributions  are 
very  smooth.  Furthermore,  the  windward  surface  has  reached  an  almost  self-similar 
behavior,  while  the  lee  side  is  still  undergoing  substantial  changes.  The  corresponding 
axial  and  circumferential  distributions  of  skin-friction  coefficient  are  shown  in  Figs.  6 
and  7,  respectively. 

The  axial  and  crossflow  distributions  of  wall  heat-transfer  rates  are  shown  in  Figs.  8 
and  9.  The  corresponding  experimental  predictions  of  the  wall  heat-transfer  rates  are 
also  shown  in  these  figures.  Figure  8  shows  the  axial  distribution  of  wall  heat-transfer 
rates,  and  shows  that  the  numerical  and  experimental  predictions  are  in  excellent  agree¬ 
ment  along  the  <f)  =  0  and  <£  =  90  deg  planes.  The  leeward  distribution  shows  a  sudden 
rise  in  the  experimental  heat-transfer  rates  along  x/Rn=6.  This  rise  is  because  the  ex¬ 
perimental  flow  on  the  lee  side  became  turbulent  around  this  region.  The  windward  side, 
however,  remained  laminar.  This  can  also  be  seen  from  the  crossflow  distribution  shown 
in  Fig.  9.  This  figure  shows  that  the  numerical  predictions  are  in  excellent  agreement 
with  the  experimental  data  up  to  (/>=  160  deg.  Beyond  this  point  the  experimental  data 
are  higher  due  to  transitional  flow  conditions. 

6.1.2.  Shock-Shape  Predictions 

The  bow  shock  shape  predicted  by  the  Case  1  calculations  is  shown  in  Figs.  10  and 
1 1.  Figure  10  shows  the  axial  variation  of  the  bow  shock  shape  along  4>  -  0,  90  and  180 
degs.  The  crossflow  distribution  of  the  bow  shock  shape  at  the  body  end  is  shown  in 
Fig.  11.  These  figures  show  that  even  under  such  a  large  angle-of-attack  condition  the 
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bow  shock  shape  predicted  by  the  new  fully  implicit,  crossflow-coupled  shock-fitting 
scheme  is  very  smooth  and  well  behaved.  This  is  despite  the  large  variations  in  the  shock 
standoff  distances  between  the  windward  and  leeward  sides.  This  clearly  reflects  on  the 
accuracy  of  the  new  implicit  shock-fitting  scheme. 

6.1.3.  Effects  of  Crossflow  Grid  Refinement 

We  also  did  some  grid-refinement  studies  to  evaluate  the  effects  of  crossflow  grid  re¬ 
finement  on  the  accuracy  and  stability  of  the  solution  scheme.  Three  different  grids  were 
used  with  31,  25  and  17  crossflow  planes,  respectively.  Figure  12  shows  the  effects  of 
this  crossflow  grid  refinement  on  the  wall  heat-transfer  rates.  This  figure  shows  that 
while  the  windward  side  remains  almost  unaffected,  the  lee  side  is  quite  sensitive  to  the 
grid  refinement.  Interestingly  the  coarse  17  plane  grid  shows  a  monotonic  decrease  of 
wall  heat-transfer  rate  from  the  windward  to  the  leeward  side.  The  25  and  3 1  plane  grids, 
however,  show  a  minimum  around  the  crossflow  separated  region,  followed  by  a  subse¬ 
quent  rise  toward  the  leeward  side.  Furthermore,  the  31  plane  solution  shows  that  the 
numerical  predictions  of  the  location  and  magnitude  of  the  minimum  wall  heat-transfer 
rate,  are  in  very  good  agreement  with  the  corresponding  experimental  data 

Figure  13  shows  the  crossflow  distributions  of  the  crossflow  skin-friction  coefficient. 
These  results  show  that  actually  the  coarse  17-plane  solution  does  not  even  predict  any 
crossflow  separation,  and  this  explains  the  corresponding  monotonic  decrease  in  the  wall 
heat-transfer  rate.  The  25-plane  and  31 -plane  results,  however,  clearly  show  a  region 
of  crossflow  separation  between  <£  =  160  deg  to  $  =  180  deg.  The  extent  of  the  separated 
region  predicted  by  these  two  fine  grids  is  also  in  good  agreement. 
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These  crossflow  grid  refinement  studies  show  that  the  grid-refinement  capabilities  of 
our  3-D  PNS  scheme  are  very  good.  Furthermore,  the  results  bring  forth  an  important 
point  that  in  order  to  maintain  acceptable  solution  accuracy  it  is  essential  to  use  ade¬ 
quate  grid- refinement  in  the  crossflow  direction.  Coarse  grids  may  produce  smooth  re¬ 
sults,  but  these  results  may  not  accurately  represent  various  flowfield  details;  such  as, 
crossflow  separated  regions,  etc. 


6.2.  Equilibrium-Air  Calculations 


The  current  equilibrium-air  calculations  were  done  in  order  to  demonstrate  the 
equilibrium-air  capability  developed  in  this  phase  I  effort.  In  terms  of  the  gas  chemistry 
effects,  the  test  case  chosen  corresponds  to  relatively  low  shock-layer  temperatures. 
Peek  shock-layer  temperatures  in  the  afterbody  region  are  of  the  order  of  400-450 
Kelvin.  Under  these  conditions  the  gas  in  the  shock  layer  behaves  like  a  perfect  gas,  and 
the  effects  of  chemical  reactions  are  negligible.  Thus,  even  if  we  use  an  equilibrium 
chemically-reacting  air  model,  if  correctly  done,  we  should  not  find  any  significant  dif¬ 
ferences  from  the  corresponding  perfect-gas  results.  Indeed  the  results  of  our 
equilibrium-air  calculations  confirm  these  observations.  Since  the  current  equilibrium- 
air  calculations  are  for  the  purpose  of  demonstrating  the  capability,  the  equilibrium-air 
calculations  were  only  done  up  to  an  axial  location  of  lORn. 

Figures  14  thru  19  show  the  effects  of  the  equilibrium-air  gas  model  on  the  surface 
measurable  quantities.  These  figures  also  show  the  corresponding  results  of  the 
perfect-gas  calculations  (Case  1).  Figures  14  and  15  show  the  axial  and  crossflow  dis- 
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tributions  of  the  wall  pressures,  while  the  axial  and  crossflow  distributions  of  the  skin- 
friction  coefficient  are  shown  in  Figs.  16  and  17,  respectively.  As  can  be  seen  from  these 
figures,  there  is  negligible  difference  between  the  perfect-gas  and  equilibrium-air  pred¬ 
ictions  for  this  case.  The  axial  and  crossflow  distributions  of  the  wall  heat-transfer  rates 
are  shown  in  Figs.  18  and  19,  respectively.  Again,  it  can  be  seen  that  the  perfect-gas  and 
equilibrium-air  predictions  for  this  test  case  are  almost  the  same.  It  should  be  noted  that 
this  close  agreement  between  the  perfect-gas  and  equilibrium-air  calculations  is  because 
of  the  low  shock- layer  temperatures  for  this  wind-tunnel  test  case.  Under  higher  tem¬ 
peratures,  more  realistic  of  free-flight  conditions,  these  differences  may  be  larger. 

Figures  20  and  21  show  the  axial  and  crossflow  variations  of  the  sock-standoff  dis¬ 
tances  predicted  by  these  perfect-gas  and  equilibrium-air  calculations.  These  shock- 
standoff  distances  are  the  radial  distances  between  the  body  and  the  shock.  Again  we 
see  negligible  differences  between  the  perfect-gas  and  equilibrium-air  calculations.  Typ¬ 
ically,  when  the  shock-layer  temperatures  are  higher,  the  equilibrium-air  effects  are 
larger  and  the  equilibrium-air  shock-standoff  distances  are  smaller  than  the  correspond¬ 
ing  perfect-gas  predictions. 

The  effects  of  the  gas  model  on  the  crossflow  separated  region  on  the  leeward  side 
are  shown  in  Figs.  22  and  25.  Figure  22  shows  the  leeward  separation  bubble  for  the 
perfect-gas  case,  while  the  separation  bubble  for  the  equilibrium-air  case  is  shown  in  Fig. 
23.  As  can  be  seen  the  location  of  separation  and  the  extent  of  the  separated  region 
predicted  using  the  two  gas  models  is  in  close  agreement.  The  detailed  pressure  contours 
on  the  lee  side  for  the  perfect-gas  calculations  are  shown  in  Fig.  24,  while  the  corre¬ 
sponding  equilibrium-air  predictions  are  shown  in  Fig.  25.  These  pressure  contours  in 
Figs.  24  and  25  correspond  to  the  same  pressure  values.  Note  the  low  pressure  region 
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around  the  separation  bubble  on  the  leeward  side.  These  figures  show  that  the  details 
of  the  crossflow  separated  region  predicted  by  the  two  gas  models  are  also  in  close  * 

agreement  with  each  other.  The  perfect-gas  calculations  predict  slightly  higher  pressures 
in  the  crossflow  separated  region;  however,  these  differences  are  still  very  small  in  mag¬ 
nitude.  ^ 

The  general  details  of  the  crossflow  solution  are  shown  in  Figs.  26  thru  28.  Figure 
26  shows  the  pressure  contours,  Fig.  27  shows  the  temperature  contours  and  the  density  * 

contours  are  shown  in  Fig.  28.  The  pressure  contours  show  the  rapid  flow  expansion 
as  we  move  from  the  windward  to  the  leeward  side.  The  temperature  and  density  con¬ 
tours  show  clearly  the  rapid  thickening  of  the  viscous  layer  in  the  crossflow  separated  ! 

region.  Again,  it  can  be  seen  that  the  flowfield  char  „ter  predicted  by  the  two  models 
is  in  close  agreement. 


6.3.  Computing  Times 

I 

The  computing  times  for  the  various  calculations  done  in  this  study  are  shown  in 
Table  2.  It  should  be  noted  that  the  computing  times  presented  in  this  table  are  for  an 

I 

IBM  3090  machine  with  scalar  LEVEL  =3  optimization.  These  results  show  that  the 
overall  computing  times  of  our  3-D  PNS  scheme  are  very  reasonable.  Since  the  solution 
scheme  involves  an  imlicit  inversion  algorithm  to  be  used  in  the  crossflow  direction,  the  ^ 

computing  times  should  increase  more  rapidly  than  the  corresponding  increase  in  the 
crossflow  grid  refinement.  Even  so,  the  increase  is  not  that  much.  For  example,  com¬ 
pared  to  Case  3,  Case  2  shows  4%  increase  in  computing  time  per  grid  point.  On  the  I 
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other  hand,  Case  1  shows  only  a  14%  increase  in  computing  time  per  grid  point.  The 
equilibrium-air  (Case  4)  computing  times  show  a  20-25%  increase  over  the  correspond¬ 
ing  perfect-gas  (Case  1)  calculations.  Although  this  increase  is  quite  reasonable  in  view 
of  the  table-look-up  procedure  used  for  the  equilibrium-air  calculations,  it  is  still  slightly 
more  than  what  we  have  typically  observed  in  our  earlier  work  (Bhutta  et  al.,  1985a). 
Typically,  we  have  observed  that  with  our  3-D  PNS  scheme  the  equilibrium-air  calcu¬ 
lations  only  take  10-15%  more  time  than  the  corresponding  perfect-gas  calculations. 
This  increase  in  the  current  equilibrium-air  computing  times  is  primarily  due  to  the  fact 
that  the  current  3-D  PNS  code  is  of  a  research  nature  and  has  not  been  fully  optimized. 
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VII.  CONCLUSIONS 


In  this  Phase  I  study  a  new  three-dimensional  and  fully-iterative  PNS  scheme  has 
been  developed  to  study  three-dimensional  hypersonic  flows  around  axisymmetric 
ballistic  configurations  under  hrge  angle-of- attack  conditions.  The  flow  around  a 
sphere-cone  configuration,  at  a  Mach  number  of  10.6  and  an  angle  of  attack  of  20  deg, 
was  predicted  to  study  the  accuracy  and  efficiency  of  this  new  three-dimensional  PNS 
scheme  under  large-angle-of-attack  conditions.  Both  perfect-gas  and  equilibrium-air  gas 
models  were  used.  The  perfect-gas  calculations  were  done  using  three  different  compu¬ 
tational  grids,  while  the  equilibrium-air  calculation  was  done  using  only  a  fine  grid. 
These  results  substantiate  the  following  conclusions: 

(1)  A  new  three-dimensional  perfect-gas/equilibrium-air  PNS  scheme  has  been  devel¬ 
oped.  This  scheme  is  an  extension  of  the  axisymmetric  perfect-gas  PNS  scheme  of 
Bhutta  and  Lewis  (1985a).  The  basic  axisymmetric  PNS  scheme  and  the  present 
three-dimensional  PNS  scheme  are  inherently  stable  in  the  subsonic  as  well  as  the 
supersonic  flow  regions  and  do  not  require  the  use  of  any  sublayer  approximation. 
Furthermore,  the  schemes  permits  very  fine  grids  to  be  used  in  the  near-wall  region 
for  improving  solution  accuracy. 


VII.  CONCLUSIONS 


42 


(2)  The  present  PNS  scheme  uses  a  second-order  accurate  smoothing  approach  which 
is  an  extension  of  the  earlier  axisymmetric  approach  of  Bhutta  and  Lewis  (1985a,b). 
In  this  approach  the  crossflow  smoothing  effects  are  applied  to  all  variables;  how¬ 
ever,  the  smoothing  effects  in  the  axis-normal  direction  are  limited  only  to  the 
pressure  field.  This  results  in  accurate  wall  heat-transfer  and  skin-friction  pred¬ 
ictions  even  with  coarse  grids  in  the  axis-normal  direction. 

(3)  A  new  predictor-corrector  solution  scheme  has  been  developed  to  treat  the  strong 
crossflow  coupling  effects  in  and  around  the  crossflow  separated  regions.  This 
predictor-corrector  scheme  involves  the  same  amount  of  operations  as  an  Approxi¬ 
mate  Factorization  scheme,  however,  it  retains  a  stronger  coupling  between  the 
crossflow  and  body-normal  directions. 

(4)  At  the  shock  a  new  fully-implicit  shock-prediction  scheme  has  been  developed  and 
used.  This  scheme  uses  a  general  curvilinear  coordinate  system  and  predicts  the 
correct  shock  location  without  having  to  make  any  approximation  about  the  viscous 
or  inviscid  nature  of  the  flow  behind  the  shock.  Furthermore  this  shock-fitting  sol¬ 
ution  is  fully  coupled  in  the  crossflow  direction  and  results  in  smoother  and  more 
accurate  shock-shapes,  and  has  very  good  stability  and  convergence  characteristics. 

(5)  It  is  shown  that  with  a  pseudo-unsteady  algorithm,  the  present  fully- iterative  three- 
dimensional  results  can  be  obtained  accurately  and  efficiently  without  any  signif¬ 
icant  computing-time  penalty.  Furthermore,  the  enhanced  solution  accuracy 
permits  much  larger  marching  steps  to  be  used  and  this  substantially  reduces  the 
total  computing  times. 
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(6)  Grid  refinement  studies  with  the  3-D  PNS  scheme  show  that  the  developed  PNS 
scheme  has  very  good  grid- refinement  characteristics.  Crossflow  grid-refinement 
results  show  that  the  present  scheme  responds  favorably  to  such  a  refinement,  re¬ 
sulting  in  substantially  more  accurate  solutions.  Furthermore,  the  results  clearly 
indicate  that  sufficient  crossflow  refinement  is  essential  to  accurately  capture  various 
flowfield  details  which  may  have  a  significant  impact  on  the  wall  heat-transfer  and 
skin-friction  predictions. 

(7)  A  three-dimensional  PNS  scheme  has  been  developed  to  accurately  predict 
supersonic/hypersonic  flowfields  around  ballistic  configurations.  The  procedure  is 
fast-convergent,  very  efficient  and  represents  a  substantial  improvement  over  exist¬ 
ing  numerical  capabilities  for  predicting  such  flows.  Furthermore,  the  scheme  is 
capable  of  treating  various  gas  models  in  a  unified  manner.  The  scheme  shows  great 
promise  for  extension  to  study  the  supersonic/hypersonic  flowfields  around  realistic 
three-dimensional  ballistic  configurations  with  fins  and  other  control  surfaces. 


l 
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APPENDIX  A.  DERIVATION  OF  3-D  PNS  EQUATIONS 


The  general  motion  of  viscous  compressible  fluids  is  described  by  the  well  known  full 


Navier-Stokes  (NS)  equations  (White,  1974).  If  we  assume  (a)  Newtonian  fluid  behav¬ 

ior,  (b)  Stokes'  Hypothesis,  and  (c)  no  body  forces,  we  can  write  the  three-dimensional 

NS  equations  as 

(pu)iX  +  (pv),  y  +  (pw)  2  =  0 

(A.l) 

(pu2  +p)iX  +  (pvu)_  y  +  (pwu)_  z  =  a, 

(A.2) 

(puv)  x  +  (pv2  +p)t  y  +  (pwv)t  z  =  a2 

(A.3) 

(puw)  x  +  (pvw),  y  +  (pw2  +p)t  z  =  a3 

(A.4) 

(pu<D),x  +  (pv<P)i  y  +  (pwd>)  z  =  a4 

(A.5) 

where 

at  =  C2PU.X  -  (2p/3)V  .  V]fX  +  [p(u>  y  +  v  x)](  y 
+  [p(u  2  +  w  x)]fZ 

(A. 6) 

a2  =  [2pv  y  -  (2p/3)V  •  V]t  y  +  [p(v  x  +  u  y)]jX 
+  [p(v  z  +  w  y)]>z 

(A.7) 

a3  =  [2 pw,  z  -  (2p/3)V  .V]z  +  |>(wtX  4-  u  z)].x 

+  [P(Wy+V  z)]>y 

(A.8) 
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and 


a4  =  V  •  (kVT)  4-  7  »(v.  t)  +  as  (A. 9) 

Here  t  is  the  stress  tensor  for  a  Newtonian  fluid,  and  is  defined  as  (White,  1974) 

*ij  =  MlXxj  +  Uj.xj  -  (2/3)<5ijukiXk]  (A.  1 1) 

where  U[  =  u,  u2  =  v,  u3  =  w,  x,  =  x,  x2  =  y  and  x3  =  z  . 

Equation  (A.l)  corresponds  to  the  conservation  of  mass,  and  Eqs.  (A.2-A.4)  corre¬ 
spond  to  the  conservation  of  momentum  in  the  x,  y  and  z  directions,  respectively. 
Equation  (A.5)  corresponds  to  the  conservation  of  energy,  and  these  equations  [Eqs. 
(A.1)-(A.5)]  are  closed  through  the  use  of  equation-of-state  for  the  gas  mixture  written 
in  the  functional  form 

fl>,T,p)  -  0  (A.  12) 

The  form  of  '<£>',  'a  5'  and  the  functional  form  of  the  equation  of  state  depends  upon 
the  gas  model  being  used.  For  a  perfect  gas  model  these  quantities  are  defined  as 


(D  -  CpT  4-  0.5V2  (A.  13a) 

a5  =  0  (A.  14a) 

fl>,T,p)  =  (yp-pT)  =  0  (A.  15a) 

For  an  equilibrium-air  gas  model  these  quantities  are  defined  as 

cD  =  h(p,T)  +  0.5V2  (A.  13b) 
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cl  5  —  0 


(A.  14b) 


Hp.T.p)  =  [p-p(p,T)]  =  0  (A.  15b) 

The  above  equations  have  been  written  in  a  nondimensional  form,  and  the 
nondimensionalization  scheme  used  is 

ui  =  u*/a^;  p  =  p*/p!0;  T 
p  =  m*/pL>5  k  =  k*/kL;  xi 

Equations  (A.1)-(A.5)  and  Eq.  (A.  12)  can  be  combined  together  and  written  in  the 
following  vectorial  form: 

el.x  +  e2.  y  +  «3,  x  =  e(gl,x  +  g2.  y  +  gj.  z)  +  P  (A- 17) 

Using  indicial  notation  we  can  write  Eq,(A.17)  as 

(ej  —£gj).Xj  =  P  (A-18) 

Or 

k!xj  =  P  (A- 19) 

Now  consider  the  general  coordinate  transformation 

<j  =  <?j(xk)  (A- 20) 

where  the  orientation  of  our  general  curvilinear  coordinate  system  is  such  that  £,  is 
measured  along  the  body,  <*2  is  measured  from  the  body  to  the  outer  bow  shock,  and  ^ 


=  T  / T^;  p  =  p  /(p^aj; 
=  x*/Rn* 


(A.  16) 


Appendix  A.  DERIVATION  OF  3-D  PNS  EQUATIONS 


47 


3  is  the  crossflow  direction.  Thus,  derivatives  in  the  transformed  space  are  related  to  the 
derivatives  in  the  physical  space  by 

(A.21) 

If  J'  represents  the  determinant  of  the  Transformation-Jacobian  for  Eq.  (A.20);  i.e., 

J  =  Detect,  <J3)/(X i.  x2.  x3)]  (A. 22) 

we  can  write  Eq.  (A.  19)  as 

(1/J)(^uk)(kk>{j)  =  (l/J)p  (A.23) 

Equation  (A.23)  can  be  further  expanded  as 

[(l/J)(^,Xk)kk]>{.  -  kk[(l/J)(£J>Xk)]>{.  =  (1/J)p  (A.22) 

Viviand  (1974)  has  shown  that  the  Jacobian  satisfies  the  identity 

Ja  =  J(Si,a),{j  (A.  25) 

where  'a'  is  an  arbitrary  quantity.  Equation  (A. 25)  can  be  used  to  obtain 

i.xk  =  J(W.?j  (A-26) 

At  the  same  time  the  chain  rule  of  differentiation  gives 

■U  =  (WJ,{j  (A-27) 
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Thus,  we  can  see  that  by  combining  Eqs.  (A.26)  and  (A. 27)  we  can  obtain  yet  another 
identity 


The  chain  rule  of  differentiation  also  gives 


(A.28) 


[(1/J)^,xk].{j  =  -  0/J)(WJ.iP  (A-29) 

Thus,  from  Eqs.  (A.28)  and  (A. 29)  we  obtain 

C(‘/J)W.{i==0  (A'30) 

Substituting  Eq.  (A.30)  in  Eq.  (A.24),  we  see  that  the  NS  equations  in  a  general 
curvilinear  coordinate  system  can  be  written  as 

[(i/j)^y,{j=(i/j)P  (A.3i) 

If  we  use  the  notation 

fj  =  (l/JKj,xkek  (A.32) 

Sj  =  (l/J)^xkgk  (A.33) 

h  =  (l/J)p  (A. 34) 

we  can  write  the  NS  equations  in  a  general  curvilinear  coordinate  system  as 

fi.5j  =  £Si.{j  +  h  (A-35) 
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In  order  to  parabolize  Eq.  (A. 35)  we  neglect  all  streamwise  diffusion  effects.  With 
this  assumption,  the  three-dimensional  parabolized  Navier-Stok.es  (PNS)  equations  in  a 
general  curvilinear  coordinate  system  become 

fj.tj  =  c(s2.£2  +S3,{j)  +  h  (A.36) 

where 

pUj 

PuU  j  4-  £j>xp 
pvUj  +  £j>yp 
pwUj  +  zp 
(hpUj 

o 

0 

mnkkU{n  +0.5(mnljUj{[i) 
mnkkv,{n  +0.5(mn2iUj>tn) 

sn  -  (p/J)  m„kkw{n  +0.5(mn3jujifn)  (A.39) 

{mnkkT,?nCp/Pr  -I-  mnkkUjUJ>{n 
(l/3)mnjkujuki{n} 

0 

where  n  =  2  and  n=  3, respectively,  and 

h-[0,  0,  0,  0,  f-flpj.p)]1"  (A. 40) 
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APPENDIX  B.  EXPRESSIONS  FOR  THE  JACOBIAN 
MATRICES 


The  three-dimensional  parabolized  Navier-Stokes  (PNS)  equations  for  a  general  gas 
flow  in  a  general  curvilinear  coordinate  system,  at  the  j  + 1  marching  step  and  at  the 
n+  l  iteration  level,  can  be  written  in  the  following  vectorial  form: 

f  jJJ ' n+1  =  cs  jjy  n+1  +  cs  jjy  n+1  +  h i+1  * n+1  (B.  1) 

Using  a  first-order  Taylor  series  expansion  around  the  previous  iteration,  we  can  write 


fj+l.n+1  j+«.n  +  A^+I,n.Aqn+I 

S  *+'• n+1  at  s  ?'• n  +  MJ2+,,n  •  Aqn+1 
S^.h+i  ^sH-i.n  +  MJ3+‘-n.Aqn+1 

hj+l.n+l^hj+I.n  +  Aji+..n4Aqn+l 


(B.2) 


where 


Aqn+’ =qj+,’n+,-qi+5,n 


(B.3) 


The  matrices  Ao,  Aj  and  M  are  called  the  jacobian  matrices,  and  have  the  following  form: 
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Aj  =  (1/J) 


-UU; 

£j.*u  +  uj 

h  yu 

Z}.zU 

0 

-VU; 

Sj,xv 

ky  V  +  Ui 

*).zV 

0 

«l.y 

— wUj 

£i.xW 

^  j.  yw 

CfciW+Uj 

0 

j.  Z 

-(DUj  ^iXd>  +  uUj  <fj>y<D  +  vUj  ^.(D  +  wUj  Uj 
0 


Mn  =  (WJ) 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Mn2l 

Mn22 

Mn24 

0 

0 

Mn3, 

Mn32 

M„33 

^n34 

0 

0 

^n41 

^n42 

Mn43 

^n44 

0 

0 

MnS1 

Mn52 

Mn53 

^nS4 

Mn55 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

o  = 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

f.pT 

0 

0 


(B- 


(B. 


(B- 


where  the  definition  of  <D  depends  upon  the  gas  model  used.  For  a  perfect-gas  model 


d)  =  CpT  +  0.5V2 


(B.9 


For  an  equilibrium-air  gas  model 
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<D  =  h(p,T)  +  0.5V2 

The  elements  of  the  viscous  jacobian  matrices  M2  and  M3  are 


Mna,  =  -  nWu /p){n  +  (l/3)mtllJ(uj/p){o 

(B.  10) 

=  (l/3)[(3mnkk  +  mnuXl/p)^] 

(B.  11) 

Mra,  =  m™nn(llP),tn 

(B.  12) 

Mn24  =  (l/3)mnl3(l /p),{n 

(B.13) 

Mnai  =  “  mnkk (v/p)>{n  +  (l/3)mn2j(Uj/p)  jn 

(B. 14) 

M„32  =  (l/3)nWl//>),{n 

(B.  15) 

Mn33  =  (l/3)C(3mnkk  +  mnjjXl/p)^] 

(B.  16) 

Mn34  =  (l/3)mn23(l/p)i{n 

(B.17) 

M„4I  =  -  mnkk(w/p)i{n  +  (l/3)mn3j(uj/p),{n 

(B.  18) 

M„42  =  (1/3)11^3,  (l/p),?n 

(B.19) 

Mn43  =  (l/3)mn32(l/p)>{n 

(B.20) 

Mn44  =  (l/3)C(3mnkk  +  mn33Xl/p),{n] 

(B21) 

-  (T/p)j{nCp/Pr  -  mnkk(V2/p)  {n  -  2mnjkuJ(uk/p)  {n/3 

(B.22) 
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Mn52  =  nw(u/p).fn  +  mnlJ(uJ/p){n/3 

(B.23) 

Mn53  -  mnkk(v/p),{n  +  mn2j(uJ/p)  {n/3 

(B.24) 

M„54  =  mnkk  (w/p),{n  +  “nS^Uj/pW3 

(B.25) 

^n55  —  mnkk(l/P),{nCp/Pr 

(B.26) 
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APPENDIX  C.  SECOND-ORDER  ACCURATE 
SMOOTHING  TERMS 


The  governing  three-dimensional  PNS  equations  are  elliptic  in  the  £2  and  <*3  directions 
so  that  we  use  central-differenced  approximations  for  all  £2  and  <*3  derivatives.  However, 
as  was  also  noted  by  Schiff  and  Steger  (1979),  the  use  of  central-differenced  schemes  is 
typically  associated  with  solution  oscillations.  This  oscillatory  behavior  becomes  more 
pronounced  if  the  local  velocities  are  small,  so  that  the  diagonal  terms  of  the  jacobian 
matrices  become  relatively  small  also.  In  order  to  damp  these  solution  oscillations,  it  is 
necessary  to  add  some  additional  higher-order  diffusion  terms  to  the  governing  PNS 
equations. 

The  search  for  an  appropriate  form  of  the  higher-order  diffusion-like  terms  which 
would  permit  a  simple  yet  a  fully  consistent  and  fully  implicit  treatment  was  very  tedious 
and  involved.  The  use  of  central-difference  formulas  for  <j;2  and  <J3  derivatives  makes  the 
solution  of  the  PNS  equations  second-order  accurate,  that  is  to  say  the  leading  trun¬ 
cation  error  is  0(A^§,  A^).  Thus,  if  we  were  to  add  O(A^)  and  O(A^)  diffusion-like 
terms  to  the  right-hand  side  of  governing  equations,  we  would  not  affect  the  formal 
second-order  accuracy  of  the  difference  scheme  in  the  £2  and  £3  directions.  The  gov¬ 
erning  equations  can  thus  be  written  as 

^sj^  +  ss^  +  h*1  +  [ni(qJ+1)](A^)2  +  Cn2(qj+1)](A^3)2  (C.l) 
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The  proper  choice  of  smoothing  terms  was  actually  based  on  a  trial  and  error  proce¬ 
dure.  To  start  with,  an  explicit  relation  relating  the  smoothed  (q)  and  unsmoothed  (x) 
variables  was  chosen  such  that  it  included  some  second-order  diffusion  effects.  A  back¬ 
tracking  approach  was  used  to  obtain  the  corresponding  smoothing  terms  that  needed 
to  be  added  to  the  governing  equations  to  produce  the  desired  result.  Once  these 
smoothing  terms  in  the  governing  equations  were  obtained,  the  order  of  each  of  these 
terms  was  analyzed  and  lower  order  terms  were  eliminated.  Then,  the  governing 
equations  with  the  modified  smoothing  terms  were  analyzed  to  see  the  impact  of  these 
changes  on  the  final  relationship  between  the  smoothed  and  unsmoothed  variables. 
After  several  iterations  of  this  trial  and  error  procedure,  we  were  able  to  find  a  proper 
choice  of  these  smoothing  terms  such  that  not  only  a  second-order  accuracy  was  re¬ 
tained,  but  a  simple  and  explicit  transformation  between  the  unsmoothed  and  smoothed 
variables  was  also  retained. 

By  choosing  the  form  of  the  vector  n  to  be 

~  ~Ao)  •  +  C(Aj  ~£M2)  •  q,  {2{23,{2  +  [(A3  — eM3)  •<!>{2{2],{3}^C-2a) 

7r2(*l)  =  {(^i/A{!  — Aq)  •  +  C(A2  — eM2)  •  A, {3{j],{2  +  [(A3  — eM3)  •  q.{3{j],{jUC.2b) 

one  can  re-write  Eq.  (C.l)  as 

[fj  +  A, . (  - q .^Ai^-q^A^)]*' 

-  *0:  +  M,  •  ( -  ^(1/4  -  q, 

+e[Sj  +  Mj  •  (  —  q  —  q  j/4]]1^' 

+  [h  +  A„  •(  -  q.iAAd  -  q.,,,  A^/4)]I+1  +  0(A{j,  Ajj) 
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Now,  let  us  define  a  new  intermediate  quantity  xi+1  as 

XJ+‘  =  q*1  _  q>  {2JjA£22/4  -  q,  {}{}A£2/4  (C.4) 

So  that 

X*’  -  qj+l  =  — q  {j{2A^/4  -  q,  {j{jA^/4  =  0(A&  M])  (C.5) 

and 


(Xi+1  -  qi+1)2  =  0(A&  A{$,  A£22A£2)  (C.6) 

Using  Eqs.  (C.5)  and  (C.6)  it  can  be  shown  that  to  second-order  accuracy  we  can  also 
write  Eq.  (C.4)  as 

,i«  _  xi«  +  x  ^  j/4  +  x  (C.7) 

Now  consider  the  Taylor  series  expansion  of  vector  f,(x1+1)  around  q)+1,  i.e., 

f2(Xi+l)  =  fa(qi+')  +  •  (Xi+1  -  qJ+1) 

+  [fj,qq]i+1  *  (Xi+1  —  qi+')2  + -  (C.8) 

-fj(qJ+')  +  A,  •  (  -  q,  {2j2A£2/4  -  q,  {}{jA£2/4) 

Thus,  to  second-order  accuracy  we  can  write  the  above  expression  as 

fj(Xi+1)  *  fj(qj+l)  +  Aj  •  (  -  q,  {2{2Ad/4  -  q.  {J{JA* j/4)  (C.9) 

Similar  expansions  can  be  obtained  for  s  and  h  so  that  to  second-order  accuracy  in  A£2 
we  can  rewrite  Eq.  (C.3)  in  terms  of  an  intermediate  solution  x,+1  as 

m1+%t  =  +  h(xi+1)  +  0(A<J2)  (C.  10) 
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The  actual  solution  that  we  seek  at  the  j  +  1  step  is  related  to  this  intermediate  solution 
by  Eq.  (C.5).  If  we  use  /  to  denote  the  grid  points  in  the  direction  (i.e.,  £  -  1,2,3... 
LMAX)  and  use  k  to  denote  the  grid  points  in  the  £3  direction  (i.e.,  k=  1,2,3...  KMAX), 
we  can  further  express  q>+1  in  terms  of  the  intermediate  solution  xi+1  as 

(*'&'  -  (x£V,/)/4  +  (xi/)/2  +  (xi-,/)/4  (C.1  la) 

4k/  -  (x‘)i/+i/4  +  (x')£//2  +  (*■)£/-, /4  (C.  1  lb) 

Thus,  we  see  that  in  order  to  introduce  a  second-order  accurate  fully  implicit 
smoothing,  we  solve  a  block-tridiagonal  system  of  equations  [Eq.  (C.  10)]  which  is  iden¬ 
tical  in  form  to  the  differenced  form  of  the  original  PNS  equations.  However,  this  sol¬ 
ution  is  just  an  intermediate  solution  (xj+1).  and  the  final  smoothed  solution  (qj+1)  can 
be  explicitly  obtained  by  using  Eq.  (C.l  1).  It  should  be  noted  that  computationally  this 
procedure  is  no  more  involved  than  the  original  (unsmoothed)  differenced  form  of  the 
PNS  equations.  Furthermore,  another  important  feature  of  the  present  second-order 
smoothing  is  that  the  additional  diffusion  terms  are  proportional  to  and  ,  so  that 
the  magnitude  of  the  aforementioned  smoothing  automatically  decreases  with  a  de¬ 
creasing  £2  and  £3  grid  size,  while  still  successfully  damping  out  the  numerical  solution 
oscillations.  Also,  despite  its  final  simple  form,  the  present  smoothing  approach  is  more 
accurate  and  performs  considerably  better  than  the  conventional  smoothing  approaches 
of  Schiff  and  Steger  (1979)  and  Shanks  et  al.  (1979). 
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APPENDIX  D.  A  MODEL  DIFFERENTIAL/ALGEBRAIC 
MARCHING  PROBLEM 


The  system  of  equations  represented  by  the  governing  PNS  equations  is  not  a  pure 
differential  system,  it  consists  of  five  partial  differential  equations  coupled  through  a 
sixth  equation  which  is  a  purely  algebraic  relation.  In  the  present  treatment  this  set  of 
governing  equations  is  referred  to  as  a  'differential/algebraic  system'.  The  most  impor¬ 
tant  view  point  to  be  presented  in  this  section  is  that  the  character  classification  of  a 
differential/algebraic  system  is  significantly  different  from  the  classical  character  classi¬ 
fication  of  purely  differential  systems.  In  other  words  a  purely  differential  system  has  a 
certain  character;  i.e.,  it  is  either  elliptic  or  time-like  or  mixed.  By  'time-like'  we  mean 
that  a  differential  system  is  either  hyperbolic  or  parabolic  or  mixed  hyperbolic-parabolic. 
However,  as  long  as  the  differential  system  is  time-like,  the  numerical  solution  can  be 
marched  in  the  time-like  direction.  On  the  other  hand,  if  the  differential  system  is  elliptic 
in  character,  marching-like  numerical  solutions  are  invalid. 

The  case  of  the  differential/algebraic  systems  is,  however,  quite  different.  For  such 
differential/algebraic  systems  the  overall  character  of  the  system  may  depend  upon  the 
way  in  which  the  problem  is  formulated.  That  is  to  say,  it  may  be  possible  to  have  a 
differential/algebraic  system  as  elliptic  or  conditionally  elliptic  if  one  formulates  the 
problem  in  one  way,  and  have  it  unconditionally  marching-like  if  one  formulates  the 
problem  in  another  way.  This  idea  is  new  and  has  given  rise  to  a  fair  amount  of  con¬ 
troversy.  Nonetheless,  it  may  be  analytically  demonstrated  on  a  model  mixed-type  sys¬ 
tem. 
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Consider  the  following  system  involving  3  unknowns,  0„  02  and  03;  i.e., 


0i,x  —  02, 0 

$2*  ~  01, y  +  202y  +  03>y=  0  (D'l) 

a20i  -03=0 


with  initial  condition  specified  at  x  =  0  and  boundary  conditions  specified  at  y=0  and 
at  y=  1.  Suppose  we  consider  the  following  initial  and  boundary  conditions: 

03(O,y)  =  a20i(O,y) 

0i(o,y)  =  02(o.y)  =  y 

03(x,O)  =  a20,(x,O) 

0liy(x,0)  =  02>y(x,O)=l 


and 


03(x,l)  =  a20,(x,l) 
0i(x,l)  =  1  +x 
02(x,l)=  1  —  (a2  +l)x 


(D.4) 


The  solution  to  Eqs.  (D.l)  for  these  boundary  conditions  [Eqs.  (D.2),  (D.3)  and 
(D.4)]  is  (Bhutta  and  Lewis,  1985d) 

0i(x,y)  =  y +  x 

02(x,y)  =  y-(a2+l)x  (D.5) 

0  3(x.y)  =  a2(y  +  x) 


The  above  model  problem  resembles  the  inviscid  limit  of  the  governing  PNS 
equations.  The  model  problem  involves  only  first-order  derivatives  in  the  two  spatial 

Appendix  D.  A  MODEL  DIFFERENTIAL/ALGEBRAIC  MARCHING  PROBLEM  60 


coordinate  directions  to  simulate  the  convective  derivatives  of  the  inviscid  limit  of  the 
PNS  equations.  The  third  equation  of  this  model  problem  is  an  algebraic  relation,  and 
is  used  to  simulate  the  role  of  the  algebraic  equation  of  state  in  the  PNS  equations.  Just 
like  the  equation  of  state  in  the  governing  PNS  equations,  the  algebraic  relation  of  the 
model  problem  appears  not  only  as  a  relation  to  be  satisfied  within  the  solution  domain, 
but  it  also  appears  in  the  initial  conditions  and  the  boundary  conditions.  The  variable 
<f> 3  of  the  model  problem  plays  a  similar  role  as  played  by  pressure  in  the  governing  PNS 
equations.  Now  consider  the  following  two  different  formulations  of  the  model  problem. 


D.l.  Formulation  I 


In  this  approach  we  can  substitute  the  third  equation  of  Eq.  (D.l)  into  the  second 
equation  and  obtain 


02, X  +  —  l)0I,y  +  2<£2,y=  0 


(D.6) 


Or  we  may  simply  write 


<Px  +  A  •  <py  =  0 


(D.7) 


where 


<P=l>lr 


(D.8) 


and  the  eigenvalues  of  A  are 
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At  =  I  +  (2  -  a2)0  5 
X2  =  1  -  (2  -  a2)0  5 


(D.9) 


It  is  shown  by  Bhutta  and  Lewis  (1985d)  that  we  can  write  Eq.  (D.7)  as 


0U  +  ^i0i.y“° 
02,x  +  ^202, y  =  0 


(D.  10) 


where 


01  =  01  +  02 
02  =  — ^101  —  ^202 


(Dll) 


Thus,  we  see  that  Eqs.  (D.10),  and  equivalently  Eqs.  (D.6),  are  time-like  if  A,  and 
/l2  are  real  (a2^2).  When  2,  and  X2  are  not  real  (a2>2),  Eqs.  (D.6)  and  (D.10)  are  elliptic 
in  nature.  In  other  words,  a  marching-like  solution  of  Eqs.  (D.6)  will  be  valid  only  if  a2 
<2.  Furthermore,  it  has  been  shown  by  Bhutta  and  Lewis  (1985d)  that  for  a2^2  the 
analytic  solution  to  Eqs.  (D.6)  can  be  found,  and  it  is  the  same  as  given  by  Eqs.  (D.5). 


D.2.  Formulation  II 


From  the  earlier  discussion  on  Formulation  I  of  the  model  problem,  we  see  that  the 
variable  'a'  of  the  model  problem  is  like  the  speed  of  sound  in  the  governing  PNS 
equations.  That  is  to  say,  in  the  classical  PNS  schemes  where  the  speed  of  sound  ap¬ 
pears  in  the  eigenvalues  through  the  pressure  terms  and  the  accompanying  equation  of 
state  (SchifT  and  Steger,  1979,  and  Shanks  et  al.,  1979),  in  the  model  problem  the  vari- 
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f 


able  'a'  appears  in  the  eigenvalues  through  the  variable  03  and  the  corresponding  alge¬ 
braic  relation.  Now,  for  this  model  problem,  if  we  can  devise  another  formulation  such 
that  the  variable  'a'  no  longer  contributes  to  the  eigenvalues  of  the  system,  it  may  pro¬ 
vide  us  with  a  key  to  attempt  a  similar  treatment  of  the  governing  PNS  equations.  Such 
a  re-formulation  of  the  model  problem  is  mathematically  possible,  and  it  will  be  called 
Formulation  II. 

In  Formulation  I,  the  overall  differential/algebraic  system  was  reduced  to  a  pure  dif¬ 
ferential  system.  In  Formulation  II  we  attempt  to  solve  directly  the  actual 
differential/algebraic  system  [Eq.  (D.l)J,  and  look  at  the  character  of  the  resulting  sys¬ 
tem.  At  first  glance  it  does  not  seem  likely  to  be  able  to  do  that.  However,  it  is  possible 
to  do  such  an  analysis  if  one  looks  at  Eq.  (D.l)  as  the  limiting  case  of  a  small- 
perturbation  problem.  Such  an  approach  is  valid  as  long  as  the  small-perturbation 
problem  being  considered  allows  us  to  take  this  limit  without  any  singular  behavior. 

For  this  purpose,  consider  the  following  problem  (were  e  2;  0) 


4>\,x  ~  o 


02,x  —  0t,y  +202,  y  +  03,y—  0 

e<^3,x=  01  ~  03 


(D.  12) 


with 


0i(O,y)  =  y 

02(o.y)  =  y  (D.i  3) 

03(O,y)  =  a2y 
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01>y(x,O)  =  1 
02.y(x,O)  =  1 

£<£3iX(x,0)  =  a20i(x,O)  -  <£3(x,0) 


(D.14) 


and 

0,(X,1)=  1  +  X 

tf>2(x,l)=  1  -(a2+l)x  (D.15) 

£</>3iX(x,  1)  =  a^^x.l)  -  <£3(x,l) 

Thus,  the  small  perturbation  problem  being  presented  has  the  correct  initial  condi¬ 
tions,  and  the  boundary  conditions  at  y  =  0  and  y  =  1  are  consistent  with  the  governing 
equations  [Eqs.  (D.12)]. 

The  complete  solution  of  this  problem  is  given  by  Bhutta  and  Lewis  (1985d).  How¬ 
ever,  briefly  speaking,  Eqs.  (D.12)  can  be  written  as 

<px  +  A  •  <py  —  (l/e)B  •  <p  =  0  (D.16) 


The  eigenvalues  of  A  are 


=  1  +J2 

A2  =  1  —  Jl  (D.17) 

;,3  =  0 

It  is  shown  by  Bhutta  and  Lewis  (1985d)  that  we  can  write  Eqs.  (D.16)  as 
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where 


0i>x  +  ^i0i,y-'*2(f7£)  =  0 

02,x  +  '*202,y  +  'il(f/e)  =  O 

*3.x-(*l-h)m=o 


(D.18) 


f=  [a2(0,  +  02)  +  (a2  -1)03]/(A,  -  A2)  (D.19) 

and 

01  =  0!  +  02  +  03 

0  2  =  —  ^101  “  ^202  (D.20) 

03  =  03 

Since  <1,  and  are  real  [see  Eqs.  (D.17)],  we  can  see  that  Eqs.  (D.18)  are  uncondi¬ 
tionally  time-like,  and  a  marching-type  numerical  solution  of  Eqs.  (D.18)  will  be  un¬ 
conditionally  valid. 

In  order  to  answer  the  question,  "How  does  the  small  perturbation  problem  of  Eqs. 
(D.12)  relate  to  the  original  problem  of  Eqs.  (D.l)?",  we  can  see  that  under  the  limiting 
condition 

e->0+  (D.21) 

Eqs.  (D.12)  reduce  to  Eqs.  (D.l),  and  the  boundary  conditions  given  by  Eqs.  (D.14)  and 
(D.15)  reduce  to  the  actual  boundary  conditions  given  by  Eqs.  (D.3)  and  (D.4).  The 
initial  conditions  are  the  same  anyway. 
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A  question  arises  —  "Is  it  valid  to  take  the  limit  of  Eqs.  (D.18)?',  or  in  other  words 
'Does  Eqs.  (D.18)  behave  singularly  because  of  the  1/e  factor?'  To  answer  to  this 
question,  the  third  equation  of  this  system  [Eqs.  (D.18)]  shows  that  for  all  e 

fife -*3.x/tf  1-^2)  (D.22) 

In  other  words,  Eq.  (D.22)  indicates  that  'f/e'  is  always  defined,  if  ^3iX  is  defined  for 

e^O. 

The  demonstration  that  is  bounded  for  e2:0  [i.e.,  Eqs.  (D.18)  is  not  singular  when 
e  -»  0  +],  comes  from  the  actual  analytic  solution  of  Eqs.  (D.18).  It  is  shown  by  Bhutta 
and  Lewis  (1985d)  that  the  analytic  solution  to  Eqs.  (D.18)  and  (D.20)  is 

<t>  ,(x,y)  =  y  +  x 

<t> 2(x,y)  =  y  -  (a2  +  l)x  (D.23) 

03(x.y)  =  a2(y  +  x)  -  a2e[l  -  exp(  -x/e)] 

where 

e  £  0  (D.24) 

It  should  be  noted  that  we  are  marching  in  the  x  direction,  so  that  our  x  is  always 
positive  and  increasing.  Thus,  we  see  that  with  e-+0+,  Eqs.  (D.22)  does  not  appear  to 
be  singular  [or,  equivalently  Eqs.  (D.22)  does  not  appear  to  be  singular]  and,  further¬ 
more,  this  solution  seems  to  be  valid  even  for  e  =  0  (Bhutta  and  Lewis,  1985d).  Also, 
we  see  that  with  e-»0+,  the  solution  to  our  hypothetical  small  perturbation  problem  ap¬ 
pears  to  correctly  approach  the  solution  to  our  actual  model  problem;  i.e., 
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(D.25) 


<£i(x,y)  =  y  +  x 
4>2(*,y)  =  y-(a2+i)x 

4>  3(x,y)->a2(y  +  x) 

The  aforementioned  mathematical  exercise  is  only  used  to  present  the  conclusion  — 
There  exist  a  class  of  'differential/ algebraic'  system-of-equations  where  it  is  possible  to 
have  a  conditionally  time-like  behavior  if  one  formulates  the  problem  in  one  way,  and 
it  is  also  possible  to  have  the  same  problem  as  unconditionally  time-like  if  one  formu¬ 
lates  the  problem  in  a  slightly  different  manner. 

The  model  problem  (Eqs.  (D.l)]  presented  herein,  bears  considerable  similarity  to  the 
inviscid  limit  of  the  governing  PNS  equations.  The  classical  treatments  of  these  PNS 
equations  (Schiff  and  Steger,  1979;  Vigneron  et  al.,  1978;  Shanks  et  al.,  1979;  etc.)  cor¬ 
respond  to  the  Formulation  I  presented  earlier,  which  was  conditionally  time-like.  The 
present  scheme,  however,  follows  the  approach  of  Formulation  II,  which  had  an  un¬ 
conditionally  time-like  character. 

D.3.  Numerical  Solution  of  the  Model  Problem 

The  key  point  to  be  understood  before  attempting  a  solution  of  the  model  problem 
is  that  this  model  problem  does  not  have  an  asymptotic  limit;  i.e.,  the  solution  grows  as 
one  marches  on.  Also,  for  this  problem  the  equations  governing  the  propagation  of  er¬ 
rors  are  the  same  as  the  governing  equations.  Consequently,  the  errors  also  convect  and 
grow  as  the  solution  marches  away  from  the  initial  conditions.  This  in  itself  is  not  a  big 
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problem  because  if  the  problem  were  only  limited  to  this,  the  solution  and  associated 
errors  grow  at  comparable  rates  and  the  relative  error  (which  is  actually  the  important 
quantity)  does  not  change  much.  However,  the  problem  with  truncation  errors  is  that 
these  errors  not  only  convect  and  grow  (as  dictated  by  the  governing  equations),  but  at 
each  step  the  cause  for  the  generation  of  numerical  errors  remains  there  and  keeps  on 
adding  to  the  existing  errors.  Thus,  the  net  effect  of  this  is  that  the  growth  of  truncation 
errors  is  faster  than  the  actual  growth  of  the  solution  itself,  and  if  one  does  not  control 
these  errors  they  can  (and  will)  become  large  enough  to  start  affecting  the  solution. 

It  is  important  to  note  that  the  solution  in  the  y-coordinate  direction  is  not  only  the 
initial  source  of  truncation  errors,  but  also  the  subsequent  reason  for  their  increased  rate 
of  growth  (because  at  each  step  it  adds  to  the  existing  error).  Thus,  it  should  be  clear 
that  the  growth  of  these  truncation  errors  is  actually  governed  by  the  number  of  times 
the  y-solution  (and  associated  truncation  errors)  are  committed.  In  other  words,  the 
errors  will  grow  more  rapidly  if  one  takes  smaller  steps  and,  thus,  commits  a  larger 
number  of  times  the  y-solution  and  the  associated  errors.  This  could  be  wrongfully  in¬ 
terpreted  as  a  step-size  (departure)  problem,  but  the  fact  remains  that  the  step-size  has 
nothing  to  do  with  the  actual  cause  of  the  problem  which  is  simply  the  inaccuracies  in 
doing  the  numerical  solution. 

A  clear  advantage  of  an  analytic  solution  of  the  problem  (if  possible)  is  that  it  does 
not  suffer  from  solution  inaccuracies  and  presents  the  true  picture.  Furthermore,  when 
using  numerical  solutions  to  either  validate  or  contradict  analytic  analysis,  it  is  essential 
that  either  the  numerical  truncation  errors  be  eliminated  or  controlled  before  one  can 
focus  on  anything  else  as  the  potential  source  of  any  problems  with  the  numerical  sol¬ 
ution. 
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In  the  following  sections  we  will  consider  two  different  approaches  to  the  numerical 
solution  of  Eqs.  (D.l).  The  first  approach  will  use  an  explicit  Lax  method  (Anderson 
et  al.,  1984),  and  will  focus  on  controlling  the  growth  of  errors  at  a  marching  step.  In 
the  second  approach  we  will  use  an  implict  bidiagonal  scheme  and  will  focus  on  the 
elimination  of  truncation  errors.  Thus,  by  controlling/eliminating  the  truncation  errors, 
we  will  obtain  numerical  solutions  which  will  truely  reflect  the  character  of  the  governing 
equations  and  not  the  effects  of  numerical  inaccuracies. 

The  explicit  Lax  method  has  been  used  for  the  c  =  0  case  only,  while  the  bidiagonal 
implicit  method  has  been  used  for  zero  and  nonzero  values  of  e.  The  main  thrust  of  our 
numerical  study  is  the  bidiagonal  implicit  scheme  because  it  is  the  most  accurate.  On 
the  other  hand,  the  explicit  Lax  method  (although  less  accurate)  is  also  interesting  be¬ 
cause  the  solution  procedure  is  identical  to  the  one  used  for  a  classical  two-dimensional 
wave  propagation  problem,  and  the  stability  constraint  obtained  is  comparable  to  the 
classical  CFL  condition. 

D.3.1.  Explicit  Lax  Method 

This  approach  is  modeled  after  the  Lax  method  for  solving  a  linear  first-order 
hyperbolic  set  of  equations  (Anderson  et  al.,  1984,  pp.  78).  In  this  method  the  x- 
derivatives  are  forward  differenced,  so  that  the  y-derivatives  become  explicitly  known 
from  the  previous  marching  step.  At  the  same  time  the  solution  values  at  the  previous 
marching  step  are  averaged,  which  is  an  0(Ay 2)  approximation.  The  idea  for  this  ap¬ 
proach  is  that  by  avoiding  the  implicit  solution  in  the  y-direction  we  substantially  reduce 
and  limit  the  generation  of  numerical  truncation  errors.  Furthermore,  in  this  case  an 


Appendix  D.  A  MODEL  DIFFERENTIAL/ALGEBRAIC  MARCHING  PROBLEM 


69 


1:  . 


I 


I 


eigenvalue  analysis  can  be  performed  on  the  error-propagation  equations  to  obtain  a 
constraint  condition  which  theoretically  limits  the  growth  of  errors.  Since  numerically 
at  each  step  we  still  cause  a  truncation  error  (however  small),  the  actual  errors  still  grow 
in  magnitude,  but  this  growth  is  extremely  small  if  we  maintain  the  stability  constraint. 


The  important  thing  to  note  about  this  method  is  that  for  a  given  step  size  in  the 
y-direction  (Ay),  there  is  only  an  upper  bound  on  the  marching  step  size  (Ax)  and  not  a 
lower  bound  which  would  have  implied  the  classical  departure-like  behavior.  As  for  the 
upper  bound  on  the  marching  step-size,  a  similar  analysis  of  the  wave  equations  also 
gives  an  upper  bound,  which  is  the  classical  CFL  condition  (Anderson  et  al.,  1984,  pp. 
78).  Thus,  this  upper  bound  on  the  marching  step-size  is  consistent  with  the  classical 
treatments  of  true  marching  problems. 


I 


i 


If  we  denote  the  x-direction  by  the  superscript  T  and  the  y-direction  by  the  super¬ 
script  the  Eqs.  (D.l)  take  the  form 


1  oo‘ 

V 

i+lj 

0i 

i,i+t 

0i 

0  10 

02 

=  (1/2) ( 

02 

+ 

02 

—a2  0  1 

<h 

0 

0 

i.j 


'  0 

-lO" 

0i 

Ax 

-1 

2  1 

02 

0 

0  0 

03 

(D.26) 


where  we  have  used  forward-differenced  approximations  for  the  x-derivatives  and  an 
averaging  of  the  solution  at  the  previous  (i-th)  marching  step.  Furthermore,  we  make 
use  of  the  fact  that  the  solution  at  the  previous  (i-th)  marching  step  is  correctly  known 
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and  satisfies  the  following  algebraic  relations  (which  come  from  the  governing  equations 
applied  at  the  i-th  step): 


ji,  j+l  2 

03  =  a  0, 

^2  ,  i,  j— I 

03  =  a  01J 

Thus,  the  Eqs.  (D.26)  reduce  to  simply  the  following 


0i 

i+lj 

0i 

i.i+1 

0i 

=  (1/2) ( 

+ 

02, 

.02. 

.02. 

-  (k/2) 


0  -f 

0i" 

i.j+1 

01 

( 

— 

(a2-l)  2  _ 

.02. 

.02. 

i.j-1 


(D.27) 


(D.28) 


and 


0^,,j  =  a20'1+1,i  (D.29) 

where  k  =  (Ax/ Ay).  The  numerical  solution  of  the  above  system  is  simple  as  all  the 
quantities  on  the  right-hand  side  of  Eqs.  (D.28)  are  known  from  the  previous  step  and 
the  values  of  0,  i+l*i  and  02  i+1>j  can  than  be  explicitly  obtained.  Once  0, i+1,)  is  known, 
the  value  of  0j1+1-i  is  obtained  from  Eq.  (D.29). 


Suppose  we  denote  the  errors  in  4>,  and  02  by  E ,  and  E  2,  respectively.  An  error 
analysis  of  Eqs.  (D.28)  can  be  performed  if  we  assume  these  errors  to  be  of  the  form 


E‘+1  =  ei+1  exp[  -iky] 
E;+1  =  e;+I  exp[  -iky] 


(DAO) 
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and  look  at  the  amplification  of  the  amplitudes  of  these  errors  (i.e.,  e ,  and  e  2)  as  we 
march  in  the  x  direction  (Anderson  et  al.,  1984,  pp.  78).  A  Fourier  analysis  of  the  error 
equations  can  now  be  performed  and  shows  that  error  amplitudes  will  not  grow  if 

|k2(c2  +  d2)  +  dk|  <  1  (D.31) 

where 

(c  +  id)  =  1  +  (2  -  a2)0'5  (D.32) 

Further  simplification  shows  that  the  constraint  of  Eq.  (D.31)  reduces  to  the  following 
constraint  on  k  =  Ax/Ay;  i.e., 

jk-kol  <  [k2+l/(c2-Fd2)fs  (D.33) 


where 


ko  =  — d/[2(c2  +  d2)]  (D.34) 

This  constraint  condition  gives  an  upper  bound  on  the  choice  of  'k'  (k  max)  or,  equiv¬ 
alently,  an  upper  bound  on  Ax  for  a  given  choice  of  Ay. 

Figures  D.l  to  D.3  show  the  results  obtained  with  the  aforementioned  explicit  algo¬ 
rithm.  These  results  are  for  a  choice  of  a2  =  10  which  corresponds  to  kmax  =  0.2.  Two 
cases  have  been  considered.  Case  1  uses  Ay  =0.1  and  Ax  =  0.001  (k  =  0.0 1 ),  and  Case  2 
is  for  a  choice  of  Ay  =  0.05  and  Ax  =  0.001  (k  =  0.02).  In  other  words,  Case  1  uses  11  grid 
points  in  the  y  direction  and  5000  marching  steps  in  the  x  direction,  while  Case  2  uses 
21  grid  points  in  the  y  direction  and  5000  marching  steps  in  the  x  direction.  Figures  D.l, 
D.2  and  D.3  show  the  </>,,  (f>2  and  <53  solution  at  y  =  0,  respectively.  For  both  cases  the 
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solution  was  carried  out  until  x=  5.0,  at  which  location  the  maximum  error  in  both  sol¬ 
utions  was  of  the  order  of  10  ~  29.  For  these  cases  the  truncation  errors  at  the  first 
marching  station  were  of  the  order  of  10  ~  31-10  -  32.  Figure  D.4  shows  the  effect  of  vio¬ 
lating  the  stability  constraint  (for  this  case)  of  kmax  =  0.2.  It  shows  solutions  for 
Ay  =  0.1  and  for  Ax  values  of  0.001,  0.005,  0.01,  0.02  and  0.1  (11  grid  points  in  the  y  di¬ 
rection  and  5000,  1000,  500,  250  and  50  steps,  respectively,  in  the  x  direction).  These 
values  correspond  to  k  values  of  0.01,  0.05,  0.1,  0.2  and  1.0,  respectively.  It  can  be  seen 
that  for  large  values  of  k  (>  0.2)  the  maximum  errors  start  to  increase  quite  rapidly. 

D.3.2.  Bidiagonal  Implicit  Method 

The  idea  behind  this  approach  is  to  eliminate  the  numerical  truncation  errors  result¬ 
ing  from  the  floating-point  roundoff  (especially  during  the  division  operations).  Multi¬ 
plication  by  whole  numbers  and  addition/ subtraction  operations  do  not  contribute  to 
the  roundoff  errors  as  much  as  the  division  operation  by  numbers  with  decimal  fraction 
parts.  Thus,  the  most  desirable  case  would  be  to  either  have  no  division  operations  or 
at  worst  have  division  by  whole  numbers  such  as  2,  4,  5,  etc.  (because  division  operations 
for  such  numbers  can  be  done  exactly  for  almost  all  cases).  One  way  of  achieveing  this 
would  be  to  (if  possible)  do  as  much  of  the  algebra  by  hand  as  possible  (such  as  matrix 
inversions,  multiplications,  etc.)  and  adjust  the  various  soiution  parameters  such  as  to 
have  either  no  divisions  or  divisions  by  numbers  such  as  2,  4,  5,  etc.  Fortunately,  in  the 
case  of  this  model  problem  this  is  possible  to  do.  In  order  to  do  so  in  the  following 
discussion  we  will  focus  on  the  case  of  a2  =3  and  for  k  =  Ax/ Ay =0.5.  The  choice  of 
these  variables  is  solely  based  on  the  algebra  simplification  they  introduce,  and  this 
makes  the  following  analysis  simple. 
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For  the  purely  differential  problem  of  Eqs.  (D.12),  we  note  that  if  we  use  one-sided 
(forward)  difference  operators  for  the  y-derivatives,  the  solution  at  a  given  marching  step 
can  be  integrated  downwards  from  the  specified  boundary  conditions  at  the  y  =  1 
boundary  toward  the  y  =  0  boundary,  where  the  solution  can  be  obtained  using  the  im¬ 
plicit  boundary  conditions  of  Eqs.  (D.14).  This  method  would  in  general  be  O(Ay)  ac¬ 
curate;  however,  it  will  be  exact  for  solutions  having  a  linear  behavior  in  y-direction. 
With  such  forward-difference  approximations  for  the  y-derivatives  and  backward- 
difference  approximations  for  the  x-derivatives,  Eqs.  (D.12)  take  the  following  form  (for 
a2=  3  and  Ax =0.5 Ay): 


2  1  0 

1  0  1 

-3 -1(1+ a) 
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-1 
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020 
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00a 
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(D.35) 


where  a  =  e/ Ax.  These  equations  are  further  simplified  by  doing  the  associated  matrix 
inversion  and  using  the  information  that 

(1  +  a)^+,’i+1  =  30l+um  +  a$3 i+I  (D.36) 


The  use  of  Eq.  (D.36)  is  correct  because  the  solution  progresses  in  the  decreasing  di¬ 
rection  (decreasing  y),  so  that  while  solving  for  the  solution  at  (i+  l,j)  grid  point  the 
solution  at  (i+l,j+l)  grid  point  has  already  been  obtained  and  satisfies  Eqs.  (D.35) 
applied  at  the  (i+l,j+l)  grid  point.  In  other  words,  Eq.  (D.36)  is  simply  the  last 
equation  of  Eqs.  (D.35)  when  applied  at  the  (i+  I,j+  1)  grid  point.  At  the  y  =  1  bound¬ 
ary,  Eq.  (D.36)  is  simply  the  differenced  form  of  the  boundary  condition  on  <f>3  [see  Eqs. 
(D.14)). 
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With  the  abovementioned  operations  the  solution  at  (i-t-  l,j)  grid  point  can  be  written 
as 

=  0i+1'J+1  +  k,(02+,’j+1  -  <t>V)  +  M03 J+1  -  «4>'i i)  (D.37) 

^+,’j  =-k3^»;+1-)+1  +  ^+,-i+1  +2 +  -k4(30i+,'i-^’i  +  ^j+1)(D.38) 

03+I,i  =  M*’,J  +  k4^J  (D.39) 

where 

k,  =  2(  1  +  a)/(2  —  a) 
k2  =  a/(2  -  a) 

k3  =  (2  —  a)/(l  +  a)  (D.40) 

k4  =*  a(l  +  a) 
k5  =  3/(1+ a) 

The  values  of  these  coefficients  (k  ,,  k2,  k3,  k*  and  k  5)  for  a  =  0,  1  and  3  (e  =  0,  Ax 
and  2Ax,  respectively)  are  given  in  Table  D.l.  As  can  be  seen  from  this  table,  for  a  =  0 
the  coefficients  are  whole  numbers  and  no  division  operations  are  involved  while  for 
a  =  1  and  3  only  division  by  2  and  4  is  involved  which  introduces  minimal  errors  (if  any). 
Solutions  were  obtained  with  different  grids  for  the  cases  of  a  =  0  and  1.  Figures  D.5  to 
D.7  show  the  solutions  for  the  a  =  0  case,  while  Figs.  D.8  to  D.ll  show  the  solutions  for 
the  a  =  1  case. 

Figures  D.5,  D.6  and  D.7  show  the  solutions  for  </>,,  <t>2  and  (/>3  at  y  =  0  for  the  case 
of  a  =  0.  Solutions  were  obtained  with  Ax  values  of  0.0025,  0.005,  0.01,  0.025,  0.05  and 
0.10  (corresponding  to  6,  11,21,  51,  101  and  201  grid  points  in  the  y-direction  and  4000, 
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2000,  1000,  400,  200  and  100  steps,  respectively,  in  the  x  direction).  In  each  case  the 
solution  was  marched  from  x  =  0  to  x  =  10,  and  the  numerical  results  showed  that  for  all 
these  case  the  solution  error  was  exactly  zero  at  all  stages  of  marching.  This  is  indeed 
to  be  expected,  because  the  exact  solutions  are  bilinear  space  functions  [Eqs.  (D.5)),  and 
for  such  bilinear  functions  one-sided  differences  in  the  x  and  y  directions  are  math¬ 
ematically  exact.  It  is  important  to  note  that  these  solutions  have  been  obtained  by 
putting  in  a  =  0  in  the  differenced  form  of  Eqs.  (D.I2).  No  singular  behavior  or  associ¬ 
ated  numerical  difficulty  was  encountered  in  doing  so,  and  the  solution  was  the  exact 
bilinear  solution  to  the  differential-algebraic  system  of  Eqs.  (D.l).  This  confirms  the 
earlier  conclusion  that  the  Eqs.  (D.l)  are  the  true  limiting  case  of  Eqs.  (D.12)  with 
e  =  0,  and  that  this  limit  is  valid  and  can  be  taken  without  any  singular  behavior. 

The  solutions  </>2  and  <£3  at  y=0  and  for  a  =  1  are  shown  in  Figs.  D.8,  D.9  and 
D.10,  respectively.  Solutions  are  presented  for  Ax  values  of  0.10,  0.025  and  0.01,  corre¬ 
sponding  to  Ay  values  of  0.2,  0.05,  and  0.02,  respectively  (corresponding  to  6x100, 
21x400  and  51x1000  grid  points,  respectively).  The  results  show  that  the  solutions  for 
</>!  and  <f>2  are  exact,  while  the  solution  for  $3  is  not  exact.  However,  it  can  be  seen  from 
Fig.  10  that  with  decreasing  Ax  the  numerical  solution  of  </>3  starts  approaching  the  exact 
solution.  This  behavior  is  explained  by  the  fact  that  the  actual  solution  of  </>3  for  c  >  0 
has  a  nonlinear  exponential  behavior  in  the  x  direction,  which  can  not  be  properly 
modeled  by  a  simple  backward-differenced  approximation  of  the  <f>3  x  term.  However, 
this  modeling  improves  as  the  marching  step  size  (Ax)  is  decreased,  and  this  is  exactly 
the  behavior  shown  by  the  numerical  results  for  <f>3.  The  analysis  presented  in  Section 
D.2  had  shown  that  the  solution  of  and  4>2  for  nonzero  values  of  a  (or  equivalently 
e)  were  composed  of  bilinear  space  functions.  The  numerical  solutions  of  these  variables 
do  give  the  exact  bilinear  space  variations  and,  thus,  confirm  the  analytic  solution  of 
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Eqs.  (D.12)  obtained  in  Section  D.3.  Figure  D.U  shows  the  percentage  error  in  the 
solution  of  4>i  for  the  three  different  step-size  distributions  that  were  tried.  This  clearly 
shows  the  improvement  in  accuracy  obtained  by  a  decreasing  sequence  of  marching  step 
sizes. 

It  is  worth  pointing  out  that  we  also  attempted  a  solution  of  the  aforementioned 
a  =  0  case  in  which  we  numerically  inverted  the  matrices  rather  than  using  the  afore¬ 
mentioned  coefficients  (k  [  to  k  5).  In  principal  there  should  have  been  no  difference  in 
these  two  treatments;  however,  such  was  not  the  case.  We  tried  the  case  with  Ax =0.01 
and  found  that  solution  suffered  from  truncation  errors  which  grew  in  magnitude  as  the 
number  of  marching  steps  increased,  and  after  a  couple  of  hundred  steps  these  errors 
significantly  degenerated  the  solution  accuracy.  Again,  the  only  difference  between  this 
calculation  and  the  similar  calculation  using  Eqs.  (D.37-D.39)  was  the  numerical  accu¬ 
racy  of  the  solution,  and  had  nothing  to  do  with  the  character  of  the  governing 
equations  which  remain  time-like. 

D.3.3.  Concluding  Remarks  on  the  Numerical  Solution 

Based  on  the  results  of  this  study  following  conclusion  have  been  drawn. 

(a)  The  results  of  this  numerical  study  confirm  the  analysis  of  Bhutta-Lewis  model 
problem.  These  results  show  that  the  solution  of  the  Bhutta-Lewis  model  problem 
can  be  properly  and  accurately  marched,  provided  one  controls  or  eliminates  the 
numerical  errors  and  maintains  solution  accuracy. 
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(b)  The  numerical  results  confirm  that  the  model  differential-algebraic  problem  is  the 
true  limiting  case  of  a  purely  differential  system  with  an  unconditional  time-like 
character.  This  limit  can  even  be  taken  numerically  without  any  singular  behavior, 
and  this  limiting  solution  is  the  correct  solution  to  the  model  differential-algebraic 
problem. 

(c)  For  marching  problems  which  do  not  have  an  asymptotic  limit,  it  is  important  to 
control  the  growth  and  generation  of  truncation  errors.  These  truncation  errors 
not  only  convect  and  grow  in  the  same  way  as  the  actual  solution,  but  in  addition 
to  this  new  truncation  errors  are  also  generated  at  each  marching  step.  Conse¬ 
quently,  the  rate  of  growth  of  these  errors  may  be  more  than  the  rate  of  growth  of 
the  actual  solution.  Thus,  after  a  sufficient  number  of  marching  steps,  these  errors 
may  become  large  enough  to  substantially  degenerate  the  solution.  In  other  words, 
in  case  of  an  inaccurate  numerical  algorithm  the  smaller  the  marching  step-size,  the 
larger  the  number  of  times  we  have  to  do  an  inaccurate  solution  in  the  y  direction 
and,  thus,  the  larger  the  number  of  times  we  cause  the  numerical  errors  to  generate. 
This  behavior  could  wrongfully  be  thought  of  as  a  problem  caused  by  the  small 
step-size;  whereas,  the  small  marching  step  has  nothing  to  do  with  the  problem  ex¬ 
cept  that  by  doing  so  the  numerical  inaccuracies  start  affecting  the  solution  much 
sooner.  The  real  solution  to  these  problems  lies  not  in  the  step-size  manipulation, 
but  in  improving  our  solution  accuracy  and  limiting  or  eliminating  associated  nu¬ 
merical  truncation  errors. 

I 
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Table  D.l.  Coefficients  of  the  Bidiagonal  Implicit  Scheme 
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Fig.  D.2.  Solution  for  d^Cx.y^O)  with  the  Lax  method  for  the  case 
of  3^=10  and  6=0 
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Fig.  D.3.  Solution  for  d>3(x,y=0)  with  the  Lax  method  for  the  case 
of  &2=10  and  t=0 
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k=0. 01  (Ay=0. 10,  Ax=0.001,  a2=10,  c=0) 
k=0 . 05  (Ay=0. 10,  Ax=0 . 005 ,  a2=10,  £=0) 
k=0 . 10  (Ay=0. 10,  Ax=0 . 010 ,  a2=10,  £=0) 
k=0 . 20  (Ay=0.10,  Ax=0.020,  a2=10,  £=0) 
k=l. 00  (Ay=0. 10,  Ax=0 . 100,  a2=10,  £=0) 


Explicit  Lax  Method 


Fig.  D.4.  Effect  of  increasing  the  ratio  k=Ax/Ay  on  the  maximum 
absolute  error  in  <£3. 
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(x,y=0) 


o  Ax=0. 0025  (Ay=0. 0050,  k=0.50,  t=0,  a2= 3) 
A  dx=0. 0050  (Ay=0.0100,  k=0.50,  e=0,  a2=3) 
+  Ax=0. 0100  (Ay=Q. 0200,  k=0.50,  e=0,  a2=3) 
x  Ax=0.0250  (Ay=0.0500,  k=0.50,  s=0,  a2=3) 
0  Ax=0. 0500  (Ay=0. 1000,  k=0.50,  c=0,  a2=3) 
^  Ax=0. 1000  (Ay=0. 2000,  k=0.50,  *=0,  a2=3) 


Fig.  D.6.  Solution  for  ^(x.y^O)  with  the  bidiagonal  implicit  method 
for  of  a2=3  and  c-0  («=0) 
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M*.y=° 


Fig.  D.10.  Solution  for  d>3(x,y=0)  with  the  bidiagonal  implicit  method 
for  of  a^=3  and  e=Ax  («=1) 
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Fig.  D.ll.  Improvement  in  the  solution  of  <£3  by  decreasing  the 
marching  step  size  (A=x)  for  a2=3  and  c=Ax 
(a=l)  and  using  the  bidiagonal  implicit  method 
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Table  1.  Freestream  conditions  for  test  case 


Quantity 

Mach  number 

10.600 

Reynolds  number 

1.32E+5 

Pressure 

(N/mJ) 

132.061 

Density 

(Kg/m1) 

9.71276 

Temperature 

(K) 

47.3377 

Velocity 

(m/sec) 

1.46E+3 

Wall  temperature 

(K) 

300.000 

Angle  of  attack 

(deg) 

20.000 

Table  2.  Comparison  of  total  computing  times  a 


Case 

x/Rn 

From-to 

*i 

Grid 
it  f» 

Time  a 
(m:s) 

Case 

1 

1.75-17.5 

28 

x  50 

x  31 

Case 

2 

1.75-17.5 

26 

x  50 

x  25 

Case 

3 

1.75-17.5 

26 

x  50 

x  17 

2:07 

Case 

4 

1.75-10.5 

16 

x  50 

x  31 

2: 10 

(a)  On  IBM  3090  with  VS-compiler  and  scalar  LEVEL=3  optimization 


Figure  I.  Coordinate  system. 
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Figure  2.  Coordinate  transformation. 
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Crossflow  distribution  of  skin-friction  coefficient 
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Figure  8.  Axial  distribution  of  wall  heat-transfer  rates. 
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Figure  14.  Effects  of  gas  model  on  the  axial  distribution  of  wall  pressures. 
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Figure  15.  Effects  of  gas  model  on  the  crossflow  distribution  of  wall  pressures. 
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Figure  16.  Effects  of  gas  model  on  the  axial  distribution  of  the  skin-friction  coef¬ 
ficient. 
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Figure  18.  Effects  of  gas  model  on  the  axial  distribution  of  the  wall  heat-transfer 
rates. 
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Figure  20.  Effects  of  gas  model  on  the  axial  variation  of  the  bow  shock  shape. 
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Figure  23.  Crossflow  separated  region  predicted  with  an  equilibrium-air  model. 
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Figure  24.  Detailed  leeward  pressure  contours  predicted  using  a  perfect-gas  model. 
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Figure  25.  Detailed  leeward  pressure  contours  predicted  using  a  equilibrium-air 
model. 
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Figure  27.  Effects  of  gas  model  on  the  crossflow  temperature  contours, 


